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ABSTRACT 


The satellite image navigation system for AVHRR (Advanced Verj’ High Resolution 
Radiometer) imagery at the Naval Postgraduate School, referred to as Avian, has been 
modified from an operator interactive procedure to an automatic navigation procedure. 
The interactive procedure, based on operator identification of discrete landmarks, has 
been replaced with a procedure which utilizes the Defense Mapping Agency’'s (DMA) 
World Vector Shoreline (WVS) as a reference. Binary shoreline images are created from 
the satellite images and correlated with the WVS reference shoreline using a sum of ab¬ 
solute differences matching technique. The correlation is performed using reference and 
search windows selected from full resolution sub-scenes of the WVS and satellite image. 

Ten images were navigated with resulting accuracies of approximately 1.3 km. The 
resultant earth location was at least as accurate as the original Avian and did not depend 
upon expertise of the operator, as the original Avian procedure does. Thus, the auto¬ 
matic Avian procedure eliminated the subjectivity inherent in the interactive landmark¬ 
ing. while reducing the amount of expertise required to perform the navigation task. 
The automatic Avian procedure does net need a landmark atlas, so image navigation can 
be done globally since the ‘.'.'orld Vector Shoreline database has worldwide coverage. 
Currently, the only subjective and interactive step remaining is the placement of the 
reference and search windows to be correlated within sub-scenes. This step can be 
completed by an operator with no p'r ' ■: experience and not effect the accuracy of the 
navigation. 
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I. INTRODUCTION 


Image registration or rectification is the process of eliminating errors or distortions 
in images so that corresponding objects in the two images correspond in size, shape and 
location. Satellite images are distorted due to Earth curvature, Earth rotation, 
spacecraft attitude, optical sensor limitations and small perturbations in satellite orbit 
(Emer\’ and Ikeda, 1984). Registration can be between two images or between an image 
and a reference source, such as a map or chart. For satellite data to be efliciently uti¬ 
lized, it must be corrected and resampled to fit a desired map projection. This procedure 
has been termed image navigation, as the image is corrected and transformed to a known 
projection. This process locates each pixel at the appropriate geographic location 
(Emery* et al., 1989). This is a fundamental requirement for the full utilization of satellite 
data, as distortions need to be removed from the images for the satellite data to be uti¬ 
lized as efficiently as possible. 

The main goal of this thesis is to develop an automatic earth location algorithm to 
be used in conjunction with the current Naval Postgraduate School's navigation system, 
referred to as Avian. Avian utilizes ephemeris data and interactive landmarking (de¬ 
scribed later in this paper) to navigate satellite imagery. The interactive landmarking 
procedure is tedious and time consuming, and its accuracy is dependent upon the level 
of training and expertise of the operator. Here the Defense Mapping Agency’s (DMA) 
digital World Vector Shoreline (WVS) is used which locates the shoreline of the world 
accurately (90''/o of all identifiable features can be located within 500 meters (Defense 
Mapping .Agency, 1988)). This reference source is compared with the shoreline observed 
in NO.A.A's (National Oceanic and Atmospheric Administration) AVHRR (Advanced 
Very High Resolution Radiometer) imagery to replace the landmarking steps currently 
used. In automating the process the World Vector Shoreline is correlated with satellite 
data, using sections of the coastline to provide a "string of control points" rather than 
just using several manually selected ground control points (GCP's). This produces more 
timely and consistently navigated imagery from a procedure which can be applied 
globally. 

Following the introduction a background on image navigation is given in Chapter 
2. The WVS data is discussed in Chapter 3. The correlation procedure and methodol¬ 
ogy developed to automaticalh navigate the imagery is outlined in Chapter 4. Chapter 
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5 and 6 discuss the experimental design and the results of the experiment, respectively. 
These chapters are followed by conclusions and appendices. Appendix A discusses the 
original Avian procedure, Appendix B gives additional details of image navigation and 
Appendix C lists the computer code used to read the WVS data. 
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II. SURVEY OF IMAGE NAVIGATION 


A. REVIEW OF REGISTRATION PROCEDURES 

Many methods of registration have been applied to satellite imageiy. Registration 
can be done by either the control point approach or the correlation approach (Davis and 
Kenue, 1978). Recent results using both methods are now reviewed. 

1. Ground Control Point Registration 

If the differences between images is any combination of translation, rotation 
and scaling the images can be registered by ground control points (GCPs) by determin¬ 
ing the positions of a minimum of two corresponding points in the images (Goshtasby 
et al., 1986). A transformation function is obtained through a least squares analysis of 
the corresponding control points. First order (linear) transformations relate the control 
points and can be used for the images if parallel lines remain parallel in the second im¬ 
age. For non-parallel and other non-linear distortions, a higher order transformation 
must be used. 

Registering images with translational differences has been examined by many 
authors (Barnea and Silverman, 19 72; Jayroe et al., 1974; Pratt, 1974; Cordan and Patz, 
1979; Leberl and Kropatsch, 1980; Eversolc and Nasburg, 1983; Goshtasby et al., 1986). 
Images with translational differences can be registered by windowing techniques which 
use a number of windows in high variance areas of one image. The corresponding win¬ 
dows are located in the second image and the center of these windows are used as GCPs 
to determine the registration parameters (Barnea and Silverman, 1972). Rotational dif¬ 
ferences can be registered by lines and segments (Stockman et al., 1982) while scaling 
differences can be handled mathematically. Labowitz and Marvin (1986) found that 
between seven and eleven GCPs were needed in image registration of Landsat imager}'. 
Orti (1981) developed optimal distributions of ground control points (at the four comers 
of the image and at areas at the right and left ed^ es approximately one quarter of the 
distance between the corners) required to obtain a given average registration error. This 
keeps the number of points to be selected to a minimum, as manual determination of 
GCPs is tedious and time consuming. 

Registration utilizing GCPs has been applied to both Landsat and AVHRR 
images. The tracking accurac} and the availability of high quality GCPs is much more 
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limited for the AVHRR (Emery et al, 1989). As this thesis pertains to AVHRR im- 
-agery, the use of GCPs with AVHRR will be examined. 

There are basically two approaches to GCP image navigation of AVHRR im¬ 
ager}’ (Emer}' and Ikeda, 1984; Emery et al., 1989). One method (the circular orbit 
method) assumes a circular spacecraft orbit based on initial ephemeris data and relies 
on known GCPs to correct for errors in Earth shape, scan geometry, satellite orbit and 
satellite attitude, to provide improvement in image to map registration. The second 
method (the ephemeris data method) relies on highly accurate satellite ephemeris data 
and only uses GCPs to correct for timing errors a satellite attitude angles. Emer}' and 
Ikeda (19S‘i) state that seven GCPs are needed with the circular orbit method to achieve 
similar accuracy as the ephemeris data method using one GCP. Both methods are ex¬ 
amined in more detail in sections to follow. 

2. Correlation Registration 

Registration can also be accomplished by correlation techniques (Jayroe et al., 
1974; Cordan and Patz, 1979; Leberl and Kropatsch, 1980; Sullivan and Martin, 1981; 
Sun, 1981; Anuta and Davallou, 1982; Crombie, 1983; Henderson et al, 1985; Anuta and 
McGillem, 1986). The two images are decomposed into a number of overlapping sub¬ 
scenes. The local shift parameters are obtained by a similarity detection algorithm after 
cross-correlation of a sub-scene from the reference image with a sub-scene from the 
other image. The shift parameters are obtained for all of the sub-scenes and a least 
squares analysis is used to compute the best fit linear transformation (Davis and Kenue, 
1978). 

Dinar}' gr:, mt images have also been utilized for image registration (Jayroe et 
al, 1974; Cordan and Patz, 1979)., These images consist of only two pixel values. For 
example, when using grey level values, the minimum and maximum grey level value may 
be the only grey level values contained at each pixel The binary images are obtained 
by thresholding the gradient images with respect to some threshold value and then the 
binar}’ gradient pictures are registered by translation, rotation and scaling of sub-scenes 
(Davis and Kenue, 1978). Anuta and McGillem (1986) discuss three correlation meth¬ 
ods of registration. These are direct correlation of gray tone image patches, correlation 
of edge images and intersecting line segments (which they termed the parameter space 
method). The parameter space method is based on the assumption that the images 
contain linear features which can be described by sets of parametrically defined line 
segments. These can be matched with the corresponding parameters of another image 
and the registration transformation determined. These three techniques were found to 
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be equivalent at high signal-to-noise ratios, but the parameter space method proved to 
be more successful at low signal-to-noise ratios. 

Edge images have been used for image navigation. Edge detection is an image 
segmentation method based on the discontinuity of gray levels or textures at the 
boundary between objects. An edge separates two regions of relatively uniform but 
different gray levels or textures. Edges are obtained by edge enhancement techniques, 
such as high pass filtering, Laplacian operators and gradient operators {.Moik, 1980). 
Nack (1977) applied edge image correlation techniques to Landsat-1 multispectral digital 
image data., Henderson et al. (1985) examined edge- and shape-guided correlation of 
control point areas and detailed way.*' of using edge descriptors in various registering 
procedures. Novak (1978) discussed correlation algorithms using edge detection for ra¬ 
dar images. Gupta and Wintz (1975) developed a boundaiy' finding algorithm for lo¬ 
cating gray level and or texture edges based on hypothesis testing. Wong and Hall 
(1979) examined several scene matching techniques, including scene matching with edge 
features which utilizes a correlation coefficient in its process. 

A specific utilization of a correlation technique is template matching, which is 
applied in the procedure de\ eloped in this paper. Template matching is the process of 
locating the position of a sub-image inside a larger image. The sub-image is called the 
template (or the window area or the reference window area) and the larger is called the 
search area (Hall, 1979; Moik, 1980; Eversole and Nasburg, 1983; Goshtasby et al., 
1984). Examples of the templates on the reference image and the corresponding search 
areas on the search image are shown in Figure 1. 

The matching process involves shifting the template over the search area and 
computing the similarity between the template and the window in the search area over 
which the template lies. Then the shift position where the largest sintilarity measure is 
obtained can be determined. Similarity measures that have been used with this tech¬ 
nique are the sum of absolute differences (Vanderburg and Rosenfeld, 1977; Hord, 1982) 
and the cross correlation coefficient (Moik, 1980; Hord, 1982; Goshtasby et al., 1984). 
Goshtasby et al. (1984) state that the sum of absolute differences is computationally fast, 
while the correlation coefficient measure is more accurate, but computationally slow. 
Hord (1982) indicates that a major advantage of the sum of absolute differences as a 
mismatch measure is that this measure grows monotonically with the number of points 
whose absolute differences have been added into the sum. This means that when looking 
for a point of best match in a given region, the process can be stopped when the sum 
of a given position has grown larger than the smallest sum obtained up to that point 
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Fig. 1. Tcniplate.s and corresponding search areas (Moik, 1980) 

from previously evaluated positions, llord also points out that an absolute threshold 
can be used in conjunction with this technique to further increase the speed of the pro¬ 
cedure. Absolute thresholding involves setting a value that is the cutoff value for the 
suinnting procedure. Then at any point, if this value is reached, the process will be 
stopped and move to the next point. 

To increase the speed of the search process for either similarity measure, a 
two-stage template technique can be used. This method first uses a subtemplate to de¬ 
termine the possible position for a match by determining the shift position which results 
in a similarity measure above some threshold value. This eliminates investing time on 
positions which show no evidence for a match. After this is accomplished, the entire 
template can be used over the good match areas to determine the best match position. 
Two-stage template matching can also be accomplished by first using a reduced resol¬ 
ution template as the subtcmplate in the first stage (Rosenfeld and Vanderburg, 1977). 
Two-stage template matching techniques have resulted in reduced computation time for 
both sum of absolute differences (Rosenfeld and Vanderburg, 1977) and for the cross 
correlation coefficient method (Goshtasby et al., 1984). 

3. Summary of image navigation 

Currently, GCPs arc being used in conjunction with ephemeris data to navigate 
satellite imagciy. The procedure and type of ephemeris data provided determine the 
number of GCPs needed to obtain accurately navigated imagery. Additional details of 
this arc provided in the following section on AVMRR image navigation. Correlation 
techniques have been used by several authors to automatically identify GCPs to be used 
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in the registration process (Nack, 1977; Ng, 1977; Davis and Kenue, 1978; Ho and Asem, 
1984). This eliminates human interaction at the GCP acquisition step, thus reducing the 
introduction of human subjectivity and the resulting variability of the navigated image 
which will occur with different operators. These automated techniques is discussed in 
greater detail in later sections. 

B. AVHRR IMAGE NAVIGATION 

As mentioned above, there are two approaches to image navigation of AVHRR 
imager}' with GCPs. These methods are referred to as the circular orbit method and the 
ephemeris data method. Both methods make use of ephemeris data and have been used 
successfully to accurately navigate satellite imager}'. Ephemeris data are the information 
that describes the location of a satellite in the geocentric reference frame. 

1. Circular Orbit Method 

The circular orbit method assumes a circular spacecraft orbit based on only 
approximate orbital parameters from initial ephemeris data provided by NOAA. Fur¬ 
ther corrections are carried out by the use of GCPs (Emer}' and Ikeda, 1984). These 
correct for errors in Earth shape, scan geometry, satellite orbit and satellite attitude. 
This method is more time consuming than the ephemeris data method as more land¬ 
marks are needed to obtain similar accuracies. This increases the amount of time spent 
with an operator selecting GCPs, unless pre-selected points are always used. There is 
also the problem of cloud cover obscuring possible landmarks. Additionally, there may 
not be much land present in the image and therefore, few landmarking possibilities. 

Legeckis and Pritchard (1976; used an algorithm with a circular orbit and a 
spherical earth, only correcting the geometric distortions due to Earth curvature. Earth 
rotation and spacecraft roll and obtained navigated imagery (NOAA-4 data) with an 
accuracy of approximately 5 kin. NOAA/NESDIS (National Environmental Satellite 
and Information Service) and C.MS (Centre Meteroligie Spaiiale) both use automated 
versions of this procedure, using a circular orbit and spherical Earth and obtain accura¬ 
cies to within 5 km (Emery et al., 1989). Emery and Ikeda (1984) used initial ephemeris 
data (assuming a circular orbit) with seven GCPs to obtain accuracies up to 1.5 km, 
while Rauste and Kuittinen (1985) navigated NOAA images with the root mean square 
error of the least squares adjustment was between 3.7 (approximately 3.4 pixels) and 6.1 
km. Ho and Asem (1986) developed a model assuming a spherical Earth and circular 
orbit, but took into account the rotation and oblateness of the Earth, following Duck 
and King (1983) and achieved accuracies as high as approximately 3 km while using only 


7 





one GCP. The single GCP was used to correct for variations in the satellite attitude and 
inclination angle. Ho and Asem (1984) also developed an automated navigation process 
by using a maximum correlation technique to identify and locate GCPs from a set of 
previously well defined GCPs. Emery et al. (1989) discusses the correction calculation 
of satellite altitude in detail, using one GCP. This correction adjusts for image dis¬ 
tortions due to the oblateness of the earth and the ellipticity of the satellite orbit. Emery 
and Ikeda (1984) used seven GCPs to make this correction, using a least-squares fit to 
allow more GCPs to be used. Scan skew' needs to be adjusted for as well, due to the fact 
that pixels are larger and do not change much away from nadir (Brush, 1988). However, 
pixels near the sub-satellite point change size rapidly. This can be corrected for by 
linearizing the scan geometr}'. Equations for linearizing the scan geometry are review'ed 
by Emer\'et al. (1989). 

2. Ephenieris Data Method 

Image navigation can be carried out utilizing an elliptical model if there is access 
to regularly updated ephemeris data. This ephemeris data is highly accurate and is 
supplied daily by U.S. Na\'y tracking stations. Several papers have been published, each 
outlining navigation models relying on this type of ephemeris data (Brush, 1985, 1988; 
Emery and Ikeda,I984; Brunei and Marsouin, 1987). Brush (1988) and Brunei and 
.Marsouin (1987) recommend a fully elliptical orbital model of relatively low eccentricity 
(Keplerian) where orbital parameters are specified by a set of accurate orbital elements 
computed from the tracking of the satellites (Emeiy et al., 1989). The orbital parameters 
needed for image navigation (Emerj’ and Ikeda, 1984) are: 

^0 the longitude of equatorial passage for the orbit of interest 

c the local value of satellite angular velocity 

H the local value of satellite altitude 

4 the equivalent equatorial passage time 

Ma the mean anomaly 

To Epoch, midnight G.MT 

oJa argument of the perigee at Epoch 

Ml the mean orbital motion 

Ml the orbital decay rate 

the longitude of equatorial passage at Epoch 
Co the eccentricity 

/ the inclination 
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The equivalent equatorial passage time (O is extrapolated using c along with v'nhemeris 
values of the mean anomaly. These orbital parameters are calculated from the 
ephemeris data, details of these calculations are given by Emery and Ikeda (1984) and 
in the summar\- article by Emery et al. (1989). 

These parameters can then be used to calculate the line and pixel locations that 
match a selected geographic map. These corrected pixel locations are calculated in x 
(position along the scan line) and y (change in scan line) to allow intercomparison be¬ 
tween corrected images. Ho and Asem (1986) call this inverse image referencing. Direct 
image referencing is the process of distorting the geographic grid to match the satellite 
image projection, thus making it difficult to compare different images. 

Emerj' et al. (1989) also discuss a method of interpolation and remapping to 
decrease the time required for the image navigation. Several resampling techniques also 
are described by Kashef and Savvchuk (19S3). Several authors (Brush, 1985, 1988; 
Brunei and Marsouin, 1987) use ephemeris methods which do not require GCPs, result¬ 
ing in accuracies ranging from 1.6 to 10.4 km. Brunei and Marsouin used clock error 
values obtained from NOA.A/NESDIS to update the time, while Brush utilized a "nudge" 
to align the coastline boundary to account for the timing error. Emery and Ikeda (1984) 
state that after their image correcting procedure there are small deviations in geography 
resulting from short-term variations in orbital parameters and fluctuations in satellite 
attitude (roll, pitch and yaw). Changes in spacecraft attitude affect the orientation of 
the radiometer, thus distorting the image. To correct for these errors Emery’ and Ikeda 
used seven GCPs using a least squares procedure, to determine how each pixel needed 
to be shifted. However, the offsets corrected by this procedure were primarily in the 
meridional direction along the path of the satellite. These meridional offsets are prima¬ 
rily due to errors in the recorded time of image data coDcction and are all of similar size 
and can be largely corrected for by using only one ground control point (Emery and 
Ikeda, 1984). With this procedure, the authors obtained maximum accuracies up to 1.5 
km. 

Emery et al. (1989) navigated images using the elliptical model to locate the 
center point of the images and then applied the circular orbit approximation to a 
spherical Earth. Then the offset along the center of the nadir track (due to the lack of 
accurate time at the receiving station) was corrected for with a "nudge" along the 
trackline to bring the coastal boundary into agreement with the image. They also navi¬ 
gated the same images with a fully elliptical model for the entire calculation, with the 
pixel locations between these two procedures having an average offset of about three 
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pixels (approximately 3 km). The authors report navigation accuracies up to 1.5 km are 
) 

possible with these techniques. 

3. Comparison of ephemeris data and circular orbit methods 

The ephemeris data method is much faster than the circular orbit method as 
there is no operator interaction required to compute the geometric corrections and 
projection resampling needed for the navigation (Emerj' et al., 1989). However, one of 
the largest sources of error is inaccurate timing at the ground station and/or on the sat¬ 
ellite. The drift error on the timing devices on TIROS-N satellites is approximately 70 
ms per day and if there is not an accurate tii \e standard to know exactly when the data 
was collected there will be errors along the satellite track. This along track shift can be 
corrected with a "nudge". This nudge moves the entire image up or down along the nadir 
trackline until the geographic features are lined up best and often a single GCP is used 
to determine this along-track nudge (Emer\‘ et al., 1989). So the ephemeris data routine 
has options for corrections to match known GCP locations to obtain information for 
this nudge. If the errors attributed to the attitude of the satellite are to be corrected, 
three GCPs are needed. Emerj’ et al. (1989) suggests this can be neglected, however, and 
still achieve accuracies on the order of one km, as the attitude of the polar orbiiers is 
controlled to better than 0.1 degrees in roll, pitch and yaw axes. 

Both of these methods can navigate images to similar accuracies, however, the 
circular orbit method needs to use more GCPs to obtain accuracies similar to the 
ephemeris data method. This increases the amount of time and human interaction 
needed to navigate the images with the circular orbit method. 

C. NFS IMAGE ^'A^TGATIOJ^• SYSTEM (AVIAN) 

The Naval Postgraduate School navigates digital imager)’ data with an ephemeris 
data system (Avian) which utilizes highly accurate ephemeris data and an interactive 
process to select landmarks. The following description of the system is taken from 
Bethke (1988). 

The digital imagery data used by this system are provided by NOAA ground receiver 
stations, such as Scripps Institution of Oceanography, and are referred to a HRPT (High 
Resolution Picture Transmission) data. These data are obtained from the Advanced 
Ver)' High Resolution Radiometer (AVHRR) sensor flown on NOAA polar orbiting 
satellites. The ephemeris data on the digital tapes represent the orbital elements of the 
satellite at OOOOZ Greenwich .Mean Time (G.MT). However, the majority of the imagery 
provided on the HRPT tapes is generated at a time T-i-OUOOZ G.MT. From the time 
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OOOOZ to T+OOOOZ, the satellite has moved west relative to the earth in its orbit (ap¬ 
proximately 15 degrees per hour). Therefore, the satellite's position, from the ephemeris 
data, must be updated to the time T, so that the satellite's position agrees with the digital 
data provided on the HRPf tape. 

The satellite's position is updated from OOOOZ GMT using the navigation model 
w’hich is designed to take the time variance of the orbital elements and the perturbative 
effects of various forces acting upon the satellite into consideration. Given the updated 
satellite position, the sub-satellite point can be determined. The satellite's 
ascending,'descending node is then calculated using the updated ephemeris data and 
known spherical geometrj' relationships. 

Once the subsatellite and the nodal points are calculated, the positions of the land¬ 
marks I dative to the geocentric Earth reference system can be found. The satellite im¬ 
ager}’ can then be "mapped" or navigated to the Earth using the reference landmarks. 
The Avian procedure is reviewed in Appendix A and the theoretical background of the 
navigation process is given in Appendix B. 

The Avian model used at the Naval Postgraduate School has provided accuracies 
between approximately 2.5 and 5 km when using NOAA AVHRR data of 1.1 km re¬ 
solution, as reported by Bethke (1988). Results of these tests and other navigated im¬ 
agery on this system indicates that the level of training of the operator will have a great 
effect on the resulting accuracies of the navigation. St. Pierre (1988) reported that in- 
creas’’^a the level of the system's operator training would probably increase accuracy of 
this s’stem to the individual pixel level. The author of this thesis has navigated imager}' 
to the pixel level, nis illustrates one of the primar}’ reasons for automating the land¬ 
marking system, as there is a need to remove the interactive, time consuming and 
somewhat subjective landmarking procedure. An automated procedure could reduce 
processing time so that the navigated images could be obtained in nearly "real time" and 
there would also be the added benefit of removing errors that are introduced into the 
system by inexperienced or less trained operators. This would remove training time for 
the landmarking procedure as well as producing more consistently navigated imagery 
an}’where on the globe. 

D. SUMMARY OF NAVIGATION ACCURACIES 

The accuracy of satellite image navigation schemes varies depending on the method 
chosen and on the implementation of that method (i.e., what number of GCPs are to 
be used, what is the distribution of the GCPs, the amount of human interaction and 
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subjectivity introduced, etc.). Table I outlines the methods discussed in this chapter and 
shows the approximate accuracy of each method. Table 1 indicates the current accura¬ 
cies obtainable and illustrates basic accuracies which the method described in this paper 
can be compared. 

Table 1. APPROXIMATE ACCURACIES OF AVHRR IMAGE NAVIGATION 


SCHEMES 


Author 

Method 

Accuracy (in km) 

Legeckis and Pritchard 1976 

Circular Orbit. 0 GCPs 

5 

Emcrv' and Ikeda 1984 

Circular Orbit, 7 GCPs 

> 1.5 

Rauste and Kuittinen 1985 

Circular Orbit, ? GCPs 

3.7 to 6.1 

Ho and Asem 1986 

Circular Orbit. 1 GCP 

3 

Brush 19S8 

Ephemeris Data, nudge 

2 to 3 

Brunei and Marsouin 1987 

Ephemeri.^ Data, 0 GCPs 

1.6 to 10.4 

Emer\- and Ikeda 1984 

Ephemeris Data, 1 GCP 

> 1.5 

Bethke 1988 

Ephemeris Data, 2-16 GCPs 

2 to 5 


Note that in Table 1, the method by Rauste and Kuitiinen does not list the number of 
GCPs used. Their method utilized GCPs but the number of GCPs used was not given 
in their paper. 
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III. WORLD VECTOR SHORELINE 


The goal of this thesis is to use the Defense Mapping Agency's World Vector 
Shoreline (WVS) as a reference base to earth locate satellite imager}'. WVS is a digital 
data file containing the shorelines, international boundaries and country names of the 
world. The absolute horizontal accuracy requirement for WVS data is that 90% of all 
identifiable shoreline features be located within 500 meters circular error of their true 
geographic positions with respect to the preferred datum (World Geodetic System). The 
shoreline data is produced from cartographic or imageiy source material. Boundaries 
and Mean High Water shorelines in the WVS data set are referenced to the World 
Geodetic System (WGS) datum and ellipsoid. 

The WVS format uses a chain-node (or chain-code or chain encoded) data structure 
to minimize redundant data storage. Each point of intersection is explicitly defined to 
eliminate gaps or overlaps between segments of adjacent features. These points are re¬ 
ferred to as nodes (endpoints of segments) and the order of points within a segment de¬ 
fines the direction, which allows identification of areas to the left and right of the 
segment (Defense Mapping Agency. 1988). 

Chain-codes arc more efficient than sequences of points that are represented by X-Y 
coordinates (Pa\lidis, 1982). Rather than digitizing an entire picture in the conventional 
wa}. ?• mesh is plai.cd over the analog picture and the vertices closest to where the curve 
crosses the lines on tlie mesh are identified. These vertices are taken to represent the 
curve and the vertices can then be encoded in an octal (eight direction) sequence by 
giving the direction from one vertex to the next (Duda and Hart, 1973). WV'S chain- 
node utilizes a common chain-code that uses eight directions, which indicate which di 
rection the feature is entering the cell of interest (Fig. 2). 

The chain-node method differs greatly from a sequential storage method. In the 
sequential storage method, individual features, as shown in Fig. 3, would be treated as 
separate pohgons, digitized independently and stored as a string of X-Y coordinates. 
The common boundar}’ between two abutting features would be digitized twice and the 
two versions may not coincide, which could create gaps and overlaps along these 
boundaries. 

In the chain-node method any common feature boundar}’ is only digitized once, so 
that the features will abut precisely and that boundar}’ will not be stored redundantly. 
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Fig. 2. Illustration of Edge codes (Defense Mapping Agency, 1988) 

The coordinates of each segment are stored with a unique alphanumeric ID so that the 
segment can be displayed separately or with other segments as a feature. 

Using \\ VS as a rcrercncc source has several advantage over other rcrcrcncc sources, 
such as charts and maps. 1 he WVS is a very accurate source of the shoreline and covers 
the entire woild. Direction of the vertices allows easy determination of the areas to the 
left and light of the segment, so that in unfamiliar areas, land and water areas are dis¬ 
tinguishable. The WVS database is extremely useful in image navigation, as it is readily 
available with worldwide coverage on a consistent geocentric datum and, as is discussed 
later, is easily presented is a binary image format. This makes it ideal for matching 
techniques involving coastlines obtained from satellite images. 
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IV. AUTOMATING AVIAN 


A. OVERVIEW OF CORRELATION PROCEDURE 

A correlation procedure has been developed and tested that utilizes the World Vec¬ 
tor Shoreline (WVS) data and the shoreline in satellite images to automatically navigate 
AVHRR images. Leberl and Kropatsch (1980) state that for automatically merging 
maps and satellite data, there is no limit to the number of possible pattern recognition 
techniques that can be used. Several authors (Ng, 1977; Davis and Kenue, 1978; Ho and 
Asem, 1984) have discussed automatic GCP acquisition. Ho and Asem's technique used 
a library of gradient values for comparison to particular landmarks. This is a disadvan¬ 
tage as there would be a great deal of time spent to develop this library and there could 
be situations where imagery could not be navigated if the known GCPs were obstructed 
by cloud cover. The papers by Ng and by Davis and Kenue both discuss correlation 
approaches to obtain GCPs. The procedure discussed by Davis and Kenue is of partic¬ 
ular interest as it is similar to the procedure used in this paper. They used binary gra¬ 
dient images in combination with a correlation technique. Windows were chosen within 
the binary images. A correlation equation would then be applied over the windows. A 
GCP would be chosen where the largest correlation value was obtained. This technique 
was used to register between two Landsat images of the same area, taken from different 
angles. 

The most commonly used matching technique is the correlation coefficient (or nor¬ 
malized cross correlation) (Hord, 1983). The disadvantage of this technique is that it is 
relatively slow, as it involves a large number of multiplications. The technique chosen 
for this experiment is the sum of absolute differences (SAD). This technique is compu¬ 
tationally more efficient than many of the cross-correlation techniques as it not only 
does not require multiplications, it also requires fewer arithmetic operations than cross 
correlation techniques. This procedure has been applied with binary shoreline images. 
The two images that are matched consist of only two gray level values. These are either 
zero (for zero gray level) or 255 (the maximum gray level). The maximum gray level 
will comprise the shoreline and the zero values the background. To apply the SAD 
procedure, each pixel value on one image is subtracted from the corresponding pixel 
value on the second image, then the absolute value is taken and summed over that entire 
image. If there is a perfect match, the SAD value is zero. By shifting the 'reference 
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window' from the reference image (WVS image) over a 'search window' from the satellite 
image, the best matching position or lowest SAD value can be found. Hord (1982) 
states that one advantage of this method is that once the sum at a given position is larger 
than the smallest sum obtained before that point, the process can be stopped for that 
position, as it will not be a best match position. This can greatly reduce the number of 
arithmetic operations performed and thus save time. 

A flow chart (Fig. 4) illustrates key steps that were used to implement the corre¬ 
lation procedure together with Avian capabilities to produce automatically navigated 
imager}’. Most of the steps described for the original Avian procedure (Appendix A) can 
still be utilized with the automatic Avian procedure. This new method eliminates the 
GCP selection procedure (Pickcm) and modifies the procedure (Twiddle) which obtains 
the timing error and the satellite attitude errors. All of the remaining steps in Avian are 
still viable. 

B. METHODOLOGY 

1. Glean and overview 

The first steps in this procedure are identical to the start of the original Avian 
navigation procedure as described in Appendix A. Glean is the process of reading all 
of the satellite image data from the AVHRR pass. Overview produces an image on a 
display monitor of the entire satellite pass at reduced resolution. This needs to be com¬ 
pleted before sub-scenes of the image can be produced. 

2. Paint a sub-scene 

After ar r-view has been created, a sub-scene of the image is 'painted'' utiliz¬ 
ing the existing procedure of Avian. The starting line and pixel number needs to 
be recorded by tne operate., as this will be needed to run the automatic matching pro¬ 
cedure. These sub-scenes can be produced for several different channels. Figures 5 and 
6 are sub-scenes of an AVHRR pass of Southern California and Baja. Through exper¬ 
imentation, it was found that a sub-scene produced using a ratio of albedo in channel 1 
to channel 2 worked well. This produced sub-scenes where the land-water interface was 
sharply defined, and tended to reduce the effect of clouds over this boundary and over 
the water just off the coast (Fig. 6). Channel 2 also worked well, producing sharp 
coastlines but clouds reflect ver}’ brightly in this channel and could cover the coastline, 
making sections of the image unusable (Fig. 5). Figure 5 shows cloud cover off the coast 
(in the channel 2 image) and Figure 6 illustrates reduced effects of cloud cover when the 
albedo ratio technique is used. 
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3. Produce a binary satellite sub-scene 

The sub-scene created by Paint undergoes several preprocessing techniques to 
prepare the image for the SAD matching. The sub-scene is thresholdcd in order to 
produce an image that contains only two gray level values, that is to produce a binary 
image. These values are zero (for the water) and 255 (for land). A value between zero 
and 255 is chosen that produces a sharp contrast at the land/water interface. In this 
study a brightness count of 45 (approxiniately n^o reflectance) has produced sharp 
coastlines for many of the images tested. All of the pixels that contain gray level values 
lower than the cutoff value are assigned a value of zero. All pixels containing gray levels 
of 45 or higlier are assigned values of 255. Tor a sub-scene created in a visible channel, 
the land areas are white and the water areas arc black. 

4. Edge enhance the sub-scene 

\\'hcn the sub-scene has been converted into a binary image the next step is to 
edge enhance the sub-scene to produce an image consisting of a coastline only (Tigs. 8 
and 9). There are many operators used for edge enhancement. The Robert's gradient 
operator (or the Robert's cross operator) was chosen and resulted with bin.aiy images 
that contained distinct coastlines. Leberl and Kropatsch (1980) recommended the 
Robert's gradient operator as it gave high peiformance with modest computing require¬ 
ments. Tills operator is rcpiescntcd mathematically in the following equation (1). 

1d+ l))' + (g(/J+ i)-g{i+ IJ))'] (1) 

where 

R(i,j) value of Robert's gradient operator at a cell or pixel 
gfi.j) a particular cell or pixel 
i line number (or row number) 

j pixel number (or column number) 


_ 

t\j + 1 

' + 1,; 

^ /+ l,;+ 1 


Fig. 7. Robert's gradient operator (Duda and Hart, 1973) 
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Figure 7 illustrates that diagonally adjacent pixels are utilized within the procedure. This 
results in a directional derivative in each direction that is approximated by simply sub¬ 
tracting adjacent elements. This step of the procedure yields a sub-scene with a binar>' 
representation of the coastline only, with a background of zero gray level values. 

By applying the thresholding and edge enhancement techniques to Figures 5 and 
6 (which are the sub-scenes created with channel 2 and the ratio of albedo of channel 1 
to channel 2, respectively) binar}’ shoreline images of both sub-scenes are created (Figs. 
8 and 9). Both sub-scenes contain areas of unobstructed shoreline w'hich could be uti¬ 
lized with correlation or matching procedures. For this study channel 2 was chosen, 
however, as it took the Paint procedure twice as long to produce a sub-scene using the 
ratio of albedo of channel 1 to channel 2 as compared to using channel 2 alone. Also, 
only small 'windows' are used to match areas on the image, so clear coastline areas can 
be chosen casih, eliminating the need to use the channel 1-2 ratio. The ratio of albedo 
method may be of more use if the window size is increased and the need to avoid clouds 
is greater. Figures 5, 6, 8 and 9 illustrate the effects of using the ratio of albedo of 
channel 1 to channel 2 and of using channel 2 alone. Figures 6 and 9 were produced 
using the channel 1-2 ratio method and Figures 5 and 8 used channel 2. When the bi- 
naiy shoreline images (Figs. 8 and 9) are examined, it can be seen that the channel 1-2 
ratio method greatly reduced the effect of light cloud cover just off of the coast. The 
channel 2 image had a greater amount of 'noise' or clutter offshore that could adversely 
effect matching results if included in window areas. This noise is produced by clouds 
which appear as non-coastal edges within the binary' image (Fig. 8). 
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Fig. 8. Edge enhanced sub-scene from ch. 2 



Fig. 9. Edge enhanced ratio of albedo sub-scene 
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5. Obtain WVS shoreline sub-scene 

The next step is obtaining a corresponding image of WVS data, wthin the same 
projection as the satellite sub-scene. The first step in generating this image is to obtain 
the set of latitude;longitude WVS points. The program to read and select the points 
from the WVS tape is given in Appendix C. The shoreline is then projected into the POS 
(Polar Orbiter Satellite) projection by utilizing the ephcmeris data that is provided on 
each HRPT tape. The latitude/longitude boundary from the satellite sub-scene is used, 
so this results in a binaiy shoreline image in the same projection and scale as the satellite 
sub-scene. An example of WVS in a POS projection is illustrated by Figure 10. 

6. Determine >vindo>vs in sub-scenes 

Selecting corresponding windows in each sub-scene is the next step in this pro¬ 
cedure. The WVS image is considered the reference image (in the literature this has also 
been termed the template image) and the satellite sub-scene is referred to as the search 
image. Each of these images will have a 'window' chosen within them where the actual 
SAD process will be calculated (Fig. 11). The reference window is smaller than the 
search window. The reference window must be smaller because the reference wndow is 
moved throughout the search window in order to determine the shift which best aligns 
the two windows. This shift is determined by how many lines and pixels the search 
window must be translated to 'match' the reference window. In order to determine what 
size the ’..’indows should be in relation to each other (to insure a match), some a priori 
information is needed. This information was obtained by performing a forward naviga¬ 
tion using no landmarks (which means that the along-track error and the satellite atti¬ 
tude angles wert ; corrected). This allows the accuracy of the images to be determined 
before any corn..-..on process took place, thus allowing the relative size of the window's 
to be determined. That is, the approximate amount of offset betw een the tw’o images 
that arc being matched needs to be known. This enables the sizes of the windows to be 
determined so that the features that are being correlated are present in each window'.. 
The average offset obtained from ten navigations was approximately 12 km. Using this 
information (Table 2, in Chapter 6), the reference window' size is chosen to be 32 X 32 
pixels and the search window to be 64 X 64 pixels. The reference window is always be 
centered within the search window. 

7. Correlate the sub-scenes 

This procedure determines the displacement betw'een the reference and search 
window's, thus resulting with the number of lines and pixels (rows and columns) the 
search window needs to be translated to align with the reference window. Figure 12 
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rig. 10. A WVS sub-scene in POS projection 

shows the entire reference and satellite sub-scenes overlaid before the along track and 
altitude angle errors have been corrected. The matching technique cliosen for this cor¬ 
relation scheme is the sum of absolute dincrcnccs (SAD). The SAD technique is given 
by: 


SAD ^ ^2^ I s[i,j) - ?•(/ -1- u.j -i- u) i (2) 

The window values are represented by s and r, with i and j representing the line (row) 
and pixel (column) numbers and u being a shift amount. The result is zero when s and 
r arc identical. 

The reference window is originally centered with the search window. As this 
matching procedure is run, the reference window is 'moved' to align within the upper left 
corner of the search window. The SAD values (line and pixel shifts) arc calculated at 
that position. The reference window is then moved one pixel at a time across the search 
window, with the SAD values being calculated at each position. After the entire 'row' 
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Fig. 11. Examples of reference and search windows 

has been completed, the reference window is moved down one line and is sequentially 
moved across that row, calculating the SAD values at each position. This continues 
until the reference window has been moved throughout the entire search window, with 
the final position being the lower right corner of the search window, The translation 
position of the reference window that has the lowest SAD values is chosen as the shift 
amount that needs to be applied to the search image to align it with the reference image. 
Figure 13 illustrates the positioning of the reference window’ within the search window 
at its first and last computational position. 

Obtaining this shift amount allows calculation of the along track error due to 
the timing error. Six lines are scanned per second with the AVHRR radiometer. Since 
this routine calculates the line and pixel shift, the line value can be used to determine the 
along track shift needed to correct for the timing error. Within the original Twiddle 
procedure, an average of the line offset of the landmarks chosen was used for this cor¬ 
rection. Using the value obtained by matching one window area is actually similar to 
averaging many GCPs, corresponding to the nuntber of pixel values for the shoreline of 
that window, since the shift obtained is an average shift for that window. 
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Fig. 13. First (a) and last (b) placement of reference nindow 

8. Iteratively correlate sub-scenes for attitude angles 

This part of the procedure determines the 'best' satellite attitude angles, based 
on new match values calculated between the reference and search windows. The value 
for the timing error has been found in the previous step. The WVS can now be re¬ 
mapped (iteratively) using the new time and varying satellite attitude angles. The atti¬ 
tude angles are first changed by 5 milliradians. There arc three angles (the roll, pitch and 
yaw) which may have positive, negative or zero (with positive being counterclockwise 
about the principal a.xis, assuming a right hand system). This .esults in a total of 27 
co.mbinations of attitude angles. For each set, the WVS has been recalculated and 
compared to the satellite binary sub-scene to obtain a SAD error. The set of attitude 
angles with the smallest error is then changed by 0.5 milliradians. The procedure is re¬ 
peated for these 27 attitude sets. The attitude angles with the smallest SAD error from 
this group is then changed by 0.1 milliradians. After these 27 comparisons are com¬ 
pleted, the attitude angles which result in the smallest error between the WVS reference 
window and the satellite search window are chosen as the attitude angles to be used to 
correct the image for distortions created by the non-zero roll, pitch and yaw. 
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9. Applying new time and attitude angles 

■ Having obtained the time and attitude corrections, they can now be applied to 
the image. This is accomplished through existing steps in Avian (Appendix A). Forw'ard 
runs a fon^’ard navigation for every 16th line and pixel value and Paint creates a sub- 
scene of the image using the updated ephemeris information. After running Forw'ard, 
the image is then mapped into one of thirteen projections using Realmap. The 
projection used depends on the needs of the user and the locations of interest. 
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V. EXPERIMENTAL DESIGN 


The automatic Avian procedure was tested on ten AVHRR passes. All ten passes 
covered the eastern North Pacific Ocean and contained portions of the North American 
continent. Before navigating automatically, the passes were navigated with the original 
Avian technique utilizing the interactive landmarking procedure. The passes were navi¬ 
gated with zero, one and then four landmarks (GCPs are referred to as landmarks in this 
procedure). All earth locations were completed by the author, which is important as 
there is an element of subjectivity within the interactive landmarking procedure. It can 
be difficult for an inexperienced operator to choose landmarks accurately and, if chosen 
inaccurately, to recognize this fact and re-select the landmark. Having the same opera¬ 
tor perform these tests may tend to reduce the variability that would be introduced by 
different operators. Accuracies will vary depending upon operator, as is sho\\'n by 
Bethke (1988). Bethke lists three major constraints influencing the selection of naviga¬ 
tion landmarks (GCPs) and the resulting accuracy of the imageiy. These deal with the 
ability of the operator to choose a 'good'' landmark, the quality of the imagery and the 
resolution of the imagery. By limiting the experiment to a single operator, the variability 
of the first constraint is minimized. 

A. METHODS OF NAVIGATION 

Tapes were navigated without landmarks, that is the timing error and the attitude 
angles were assui : to be zero. One reason for navigating the imagery this way was to 
obtain appro.\imv.ie ofl'sets before corrections were made to the time and attitude angles 
so that the size of the reference and search windows could be determined. It was nec¬ 
essary to guarantee that shoreline in the search window would also appear in the smaller 
reference window. Secondly, this gives reference earth location accuracy values to serve 
as controls for the automatic method. 

A\ian was run again on the passes using one landmark. This was done so that im¬ 
ages navigated with one landmark could be compared to imagery navigated with one 
window correlation through the automatic Avian navigation procedure. One landmark 
is often used in other studies to correct for timing errors (Table 1). 

The ten passes were then navigated using four landmarks. This approach was cho¬ 
sen since in previous papers the accuracy often improved when using more than one 
GCP. four was the optimum number of GCPs for all of the passes. More landmarks 
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were not used as some images had a good deal of cloud cover which limited the number 
of GCPs. The landmarks were chosen as carefully as possible. Twiddle displays ap¬ 
proximate error vectors when the Twiddle procedure is initiated (Appendix A). This al¬ 
lows the operator to recognize poorly selected landmarks and to return to the earlier 
Pickem procedure and re-select the landmark. However, when choosing more than one 
landmark it is much harder to obtain all of the landmarks as accurately as when a single 
landmark has been chosen. When one landmark is selected, it can be repeatedly selected 
until the error vector is small. When more landmarks are be selected concurrently, it is 
much more difficult for the operator to choose which landmarks are are causing the 
greatest effect on the error vectors. Thus, it is harder to select each of the multiple 
landmarks as accurately as a single landmark. 

The ten passes were then navigated using the automatic Avian procedure. The 
windows were selected in clear areas of the coast, as close to the centerline of the image 
as possible. However, several of the passes were either too cloudy over the center of the 
image to allow this or the land area was located to the right half of the image. Each 
navigation used one window area to obtain the corrections to the tiu;e (yielding an along 
track error) and to the roll, pitch and yaw. 

B. ACCURACY DETERMINATION OF NAVIGATIONS 

The accuracies of the individual navigations were determined by utilizing a feature 
contained in the Paint procedure of Avian. Once the updated time and attitude angles 
were obtained (cither by Pickem and Twiddle in the original Avian or by the automatic 
Twiddle in automatic Avian) sub-scenes could be created by the Paint procedure. The 
latitude, longitude of particular points (known landmark locations were used) can then 
be obtained simpl\ by moving a cursor to the point in question. This latitude/longitude 
position is determined by the mapping procedures described in Appendix B which utilize 
the updated ephemeris data. The difference between the navigated position and the ac¬ 
tual position was then determined. This difference was obtained using the following 
equation from Bowditch (19S4j for the great circle difference: 

D = cos”' [_{sinL^ x siiiLi) -f {cosLi x cosL^ x cosDLq)'] (3) 

where 

i, latitude of navigated oint 

Li true latitude of the point 

DLf, difference between longitudes 
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D distance (in degrees) 

There are more accurate methods of determining the distance between two points; 
however, as the resolution of AVHRR imager>- is 1.1 km, this method is more than ad¬ 
equate for the needs of this experiment. 

The points chosen to validate the accuracy were selected as evenly distributed over 
the pass as possible. In cases were the entire coastline was not contaminated by clouds, 
points were selected over the entire range of the coastline and were evenly spaced. In 
images with partial obstruction due to clouds, the clear coastline areas were used, again 
with even distributions of the points throughout these areas. As different numbers of 
points were used for different navigations and different images, the standard error of the 
mean was determined for each navigation and then for each method (i.e., for zero, one 
and four la’^dmarks and for automatic Avian) to statistically weight the determinations, 
as described in the data analysis section. 

Locating the points on the image to validate the accuracy of the navigation is a 
subjective step in a sense. Locating each landmark to the exact pixel can be ver>’ diffi¬ 
cult, particularly towards the lateral sections of the image. Distortions due to the sat¬ 
ellite projection and the changing areal content of the pixels (Fig. 14) make this difficult 
at times. For this reason, points were selected as evenly distribu‘’^d as possible and then 
averages of these points obtained to validate the earth locations. 
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Fig. 14. Illustralion of the effect of scan geometry on pixel size (Betiike, 1988) 









VI. RESULTS 


The results of navigating ten different AVHRR images by the original Avian pro¬ 
cedure and the automatic Avian procedure are presented in Table 2. The columns la¬ 
belled zero landmarks, one landmark and four landmarks were navigated with the 
original Avian procedure and the last column is the results of the automatic Avian pro¬ 
cedure. The date of the satellite pass is given in the first column. 

The average accuracy with zero landmarks was approximately 12 km, with the low¬ 
est accuracy of a pass being 10.35 km and the highest at 15.03 km. These values are not 
as accurate as the procedure developed by Brunei and Marsouin (Table 1) which does 
not use any GCPs. Their method does however, make use of time errors obtained from 
XOAA.'NESDIS which would eliminate much of the along track error resultant from 
clock offset. The original Avian procedure with zero landmarks does not correct for the 
time offset, thus the accuracies are not as small. The third column presents the accura¬ 
cies of the orignal Avian technique using one landmark. The average accuracy was ap¬ 
proximately 2 km, with 1.05 being the minimum and 3.17 the maximum. These values 
compare favorably with techniques developed by Ho and Asem (1984) and Emery and 
Ikeda (1984) which also utilized one GCP. The previous studies resulted with accuracies 
ranging from 1.5 to 2.0 km. Avian with one landmark produces accuracies within this 
range. The original Avian procedure utilizing four landmarks is presented in the fourth 
column. The average accuracy improved to 1.5 km, with a minimum accuracy of 1.00 
km and a maximum of 1.97 km. The last column illustrates the accuracies produced with 
the automatic Avian procedure. The average accuracy was approximately 1.3 km, with 
a minimum accuracy of 1.00 km and a maximum of 1.69 km. Automatic Avian com¬ 
pares veiy- well with the method developed by Brush (1988). Brush's method utilized a 
'nudge' to bring the coastline within an image into alignment with a reference coastline. 
This produced accuracies of 2 to 3 km. The automatic Avian procedure utilizes shoreline 
matching to update the ephemeris data and produces imagerv’ navigated close to the re¬ 
solution of the data. 
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Table 2. ACCURACIES OF AVIAN METHODS (IN KM). The columns labelled 


0, 1 and 4 landmarks were navigated with original Avian. The last column 
used the automatic Avian procedure._ _ 


Tape date 

0 landmarks 

1 landmark 

4 landmarks 

auto Avian 

05-15-86 

12.44 

2.12 

1.63 

1.00 

05-16-86 

11.45 

2.67 

1.97 

1.69 

06-17-87 

14.35 

1.96 

1.02 

1.48 

06-28-87 

12.65 

3.17 

1.00 

1.64 

07-06-87 

12.39 

1.64 

1.13 

1.37 

07-07-87 

11.57 

1.05 

1.78 

1.43 

07-14-87 

12.08 

1.13 

1.40 

1.04 

07-16-87 

15.03 

1.52 

1.45 

1.23 

07-18-87 

10.35 

1.84 

1.74 

1.13 

04-23-88 

10.91 

2.13 

1.47 

1.22 

Mean 

12.32 

1.92 

1.46 

1.32 


In order to test whether the automatic Avian procedure is more accurate than the 
original Avian with the interactive landmarking, the data sets were statistically com¬ 
pared. As stated earlier, since varving numbers of points were used to determine the 
accuracies of different cases, the standard error of the mean was determined for each 
method (Table 3). If two populations have the same variance the Student's t statistic 
can be used to test whether the true means of two populations are the same (Mi = H:) 
However, if it is not known that the variances are equal, then the problem can be ap¬ 
proximated by the Student's t, using the procedure described as the Fisher-Behrens 
problem (Hamilton, 1964). With this procedure the degrees of freedom for the Studeijt's 
t is given in equation 4. 



■ , (^ 2 ) ' 

CM 

(«]) («2) 

- « 

, (^2) 

(/7iX(«i-l)) («2X(«2-1))_ 


(4) 


Using this test, the automatic Avian procedure is not more accurate than the original 
Avian with four landmarks at 95'’/o confidence. The automatic Avian procedure is more 
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accurate than the original Avian procedure utilizing one landmark at 95% confidence 
and of course, was significantly better that the original Avian with zero landmarks. 


Table 3. MEANS AND STANDARD ERROR OF THE MEANS OF IMAGE 


NAVIGATION 


statistic 

0 landmarks 

1 landmark 

4 landmarks 

auto Avian 

•Mean 

12.32 km 

1.92 km 

1.46 km 

1.32 km 

S.E. of Mean 

0.786 km 

0.340 km 

0.262 km 

0.246 km 


Although the original Avian procedure with four landmarks produced earth lo¬ 
cations of approximately the same accuracy as the automatic Avian procedure in this 
experiment, this is somewhat misleading. As previously stated, the original Avian pro¬ 
cedure can produce accuracies that vary greatly depending upon the level of training and 
the expertise of the operator. .Many earth location studies have been completed w'here 
the accuracy varied from 2 to 5 km. The automatic Avian procedure has eliminated the 
subjectivity inherent with operator interactive acquisition of landmarks. Presently, dur¬ 
ing the navigation procedure, the interactive procedure is limited to the operator select¬ 
ing a clear portion of the sub-scene that contains shoreline to locate the reference and 
search w'indows. The timing error and attitude angle errors are then calculated auto¬ 
matically. Figures 15 and 16 shows the satellite and WVS sub-scenes overlain before and 
after the timing error and attitude errors have been determined. The WVS has been 
calculated utilizinc ■' e ephemeris data in Figure 15 without the corrections made, and 
in Figure 16 the tuiung and attitude angle errors have been corrected, illustrating that 
these corrections do improve the ephemeris data for this satellite pass. 

The results of the earth location with automatic Avian compare favorably with re¬ 
cent studies (Table I). Accuracies are approaching the resolution level of the imagery 
itself and are obtained with little human interaction. Results of the navigation with one 
landmark using original Avian provide accuracies of approximately 2 km. This confirms 
that the conclusions reached by Emery et al. (1989) that navigation can be completed 
accurately using one GCP while ignoring attitude angle errors. Figure 15 illustrates that 
the along-track error is the largest offset between the tw'o shorelines. Additional land¬ 
marks used with the original Avian procedure result in slightly more accurate navigation, 
as a result of obtaining more accurate roll, pitch and yaw errors and applying these 
corrections with the navigation procedure. The automatic Avian procedure produces 
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accurately navigated imagery while eliminating the subjective nature of the original 
Avian procedure. Accuracies equivalent or better than recent studies has been accom¬ 
plished without relying on stored sets of GCPs, as Ho and Asem (1984) or updates of 
timing errors of the satellites from outside sources, as Brunei and Marsouin (1987). 
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Fig. 15. WA’S and satellite sub-scene overlaid before correctioii.s: The WVS is de¬ 
picted as the white shoreline. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

An automatic method of navigating AVHRR imagery has been developed and tested 
successfully in this paper. This automated procedure utilizes DMA's database as 
a reference base, binarj' shoreline images and the Sum of Absolute Differences (SAD) 
matching technique to solve for errors in timing from the satellite's clock and in the 
satellite attitude angles. When these errors have been calculated the corrections are 
added to the ephemeris data and then the existing capabilities of Avian produce navi¬ 
gated imager}' through a series of coordinate transformations, resulting in imager}' with 
standard geographic coordinates. It has been shown that this automated method can 
produce navigated imager}’ at least as accurate as the original Avian procedure with four 
landmarks. The automatic method also eliminates the subjectivity inherent with inter¬ 
active landmark acquisition, so that the navigation can be performed by inexperienced 
operators and still result with highly accurate results. 

Another advantage of this automatic navigation method is that a set or library’ of 
GCPs do not need to be accumulated and stored to be used in the navigation process. 
The original Avian procedure had a set of landmarks that had to be accumulated over 
a period of time. The automatic procedure developed by Ho and Asem (1984) had to 
rely on a set of stored landmarks, based on gray level gradients. The automatic Avian 
procedure matches the shoreline from the satellite image with the shoreline presented in 
the WVS database. Because WVS covers the entire world, there is no problem in ob¬ 
taining shoreline data. 

Currently, the automatic Avian procedure takes approximately the same amount of 
time (approximately one hour) to complete as the original interactive procedure al¬ 
though the time spent interactively selecting landmarks can vary greatly. The image 
navigation in this study vi'ere performed on a DEC Microvax II and time constraints 
will var}' depending on computer system used. Selecting four or five landmarks accu¬ 
rately can be very difficult, and an operator may have to re-select landmarks several 
times to be able to produce accuracies of navigated imagery similar to that produced in 
this paper. There is also the benefit of less interactive steps with the automatic proce¬ 
dure. Once the operator has selected the window area, the automatic Avian technique 
has no other interactive steps. This allows the operator to complete other work while 
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the corrections to the time and the attitude angles are calculated. The original Avian 
’technique involved the operator interactively throughout the procedure. At the present 
time, the automatic procedure has not significantly reduced the amount of time spent 
navigating imagery, but it has reduced the amount of time that an operator is interac¬ 
tively involved in the process. 

B. RECOMMENDATIONS 

There are several procedural steps that can be improved to produce more consistent 
and timely earth location. 

1. Placement of the matching >vindo» s 

A comprehensive study to determine whether or not the placement of the refer¬ 
ence and search windows within the satellite image has an effect on the accuracy of the 
navigation w'ould extend this research. It is likely that placement of the w’indow' along 
the lateral extremes of the image where there is more distortion will produce reduced 
accuracies. Individual passes would need to be navigated several times w'ith the w’indows 
located at vaiw ing locations along the coastline to explore the impact of window location 
on accuracy. 

2. Development of multi-window techniques 

The original Avian procedure twiddles up to 16 GCPs (or landmarks) at one 
time to obtain the optimum attitude angle corrections. Development of an automated 
procedure utilizing more that one set of reference and search windows may improve the 
accuracy of the navigation over the entire pass. If the study of window placement de¬ 
termines that lateral p.icemcnt of the window increases accuracy, a multi-window tech¬ 
nique would allow placement of windows at varying distances from the centerhne, 
possibly improving accuracies. 

3. Fully automate the Avian procedure 

Development of a fully automatic navigation procedure that would produce 
accurate earth location would be very beneficial. Reduction in the amount of time that 
the procedure needs operator attention is desirable. A possible approach would be to 
map the WVS data into the POS projection using the ephemeris data provided on the 
tape. Through this work, it is known that the satellite shoreline w'ould be within ap¬ 
proximately 12 km of the shoreline from the WVS. Depending upon the number of 
w’indow areas desirable and the number of points in the WVS for that image, an arbi¬ 
trary number can be selected that would pick a point out of the WVS. This point would 
be used as the midpoint for the reference and search windows. A sub-scene would first 
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need to be created at full resolution which contained these areas. It could be thresholded 
and edge enhanced automatically. The SAD matching values would then be calculated 
between the reference and search windows. A critical valr.e would be set (this would 
need to be determined ahead of time through experimentation) and if the error obtained 
from the SAD procedure was below this critical value, this window area could be used. 
If the window area is rejected, the procedure would skip down the recommended number 
of points within the WVS data and repeat the procedure. This could be repeated until 
a satisfactory window area is found or if a multi-vindow technique is being used, until 
the number of windows specified is found. A fully automatic Avian procedure could 
produce accurately navigated imageiy- with a minimal amount of operator time involved. 
Thus, earth located satellite imageiy would be readily available for scientific use. As the 
use of digital polar orbiter imageiy increases, automated navigation techniques, like the 
method described here, will become even more critical for the effective use of NOAA and 
DMSP data. 

4 . Compare with new navigation technique 

The Department of Oceanography at the Naval Postgraduate School has re¬ 
cently acquired a satellite image navigation software package from the University of 
Miami. This navigation system was not online at the time this thesis was printed. This 
navigation system also uses coastal boundaries to navigate satellite images. A detailed 
comparison between the automatic Avian procedure and this new package needs to be 
performed in order to determine the uses of these procedures at the Naval Postgraduate 
School. .Accuracy of navigation (including the accuracy of the shoreline reference base 
used), amount of t ' spent performing navigation and user friendliness of the two sys¬ 
tems need to be co;..pared. 
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APPENDIX A. ORIGINAL AVIAN PROCEDURE 


Avian is a software package that has been implemented at NPS for navigation, cal¬ 
ibration and image analysis of polar orbiter AVHRR data. Avian is presented here in 
the same sequence as to the user at the terminal, a step by step procedure that follows 
a menu given in the IDEA (Interactive Digital Environmental Analysis) lab. The the¬ 
oretical background of the navigation process is reviewed in Appendix B. There are 
eight steps, some of which may or may not be needed during a particular navigation. 

1. Glean 

This is the first step and is required for all tapes that are to be navigated. 
Gleaning retrieves all of the satellite image data from a magnetic tape. Glean produces 
three files, an image information file (INF), a raw’ calibration data file (CAL) and an 
image data file (SCN). After Glean is run , a routine called Calib is automatically run 
that produces a file containing actual calibration coefficients (i.e., temperature coeffi¬ 
cients for channels 4 and 5 in Kelvin, albedo coefficients for channels 1 and 2 and radi¬ 
ance coefficients for channel 3). 

2. Oveniew 

This routine allows the user to view the entire image on a display monitor, at 
any channel (although at low resolution). This step is required if an image is to be 
navigated as sub-images at full resolution are created from this in the Pickem and Paint 
steps of Avian. An overview (OVR) image file is created. 

3. P’ckem 

This routine is used to supply a set of points w'hich w’ill be used to create a 
mapping from the image coordinate system to the geographic coordinate system. The 
INF, SCN and OVR files are required to run Pickem. The over\'iew' image is presented 
and the operator then chooses small areas at which known landmarks exist (points, 
headlands and sharp edges along the coast). These areas are expanded by a factor of 
four and the operator then picks the landmark location by cursor. When selecting these 
landmarks the operator must choose landmarks that have their positions stored in a li- 
brarj’ of GCPs which is stored w'itliin the program. A list of these is available for oper¬ 
ator reference. The landmarks are placed in a landmark (LND) file, to be used in the 
routine following Pickem. 
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4. Twiddle 

Twiddle is an iterative procedure which determines four unknown quantities, 
the timing error (which accounts for an along-track error) and the satellite attitude an¬ 
gles (the roll, pitch and yaw). The Twiddle procedure is showm in Figure 17. This step 
uses the INF and LND files that w'ere produced in the earlier steps. Twiddle uses posi¬ 
tions of landmarks (GCP's) which were manually selected in the Pickem procedure. The 
pixel and line number (or row and column) were saved for each selected landmark. 
Twiddle has a library of GCP's in which the latitude and longitude of each GCP is 
stored. For each GCP chosen, the latitude and longitude of the corresponding stored 
'correct' positron is converted to line and pixel numbers. This conversion utilizes the 
ephemeris data and results is a 'navigated' GCP position, which the picked landmark 
will be compared. The conversion from geodetic coordinates to satellite image coordi¬ 
nates is a reverse navigation and is detailed in section E of Appendix B. The difference 
betw’een the line number of each landmark and its corresponding navigated GCP posi¬ 
tion is then obtained. The average of all of these offsets is then calculated. This average 
line offset yields the timing offset as the radiometer on the NOAA satellites scans six 
lines per second. The offset in time is then added to the time given in the ephemeris data 
to update and time and thus correct the along-track offset. 

To obtain the corrections for the satellite attitude angles several steps are taken. 
The first step is to calculate the error for each landmark (i.e., the difference in the posi¬ 
tion of the landmark GCP and the navigated GCP). These errors are calculated and an 
error vector is saved and also displayed on the monitor for each GCP. At this point, 
any landmark which appears to ha\'c a large error can be erased (using the procedure 
described under Utilities) and re-selected in the Pickem procedure. These error vectors 
are first obtained assuming zero errors in the roll, pitch and yaw. The landmarks are 
then 'twiddled', that is each of the three attitude angles are changed by fixed amounts 
in both directions. There are a total of 27 combinations of differing roll, pitch and yaw 
that are calculated (including the position of zero error in the roll, pitch and yaw). The 
'best' roll, pitch and yaw combination is chosen by selecting the attitude angles w'hich 
produce the minimum offset between the landmarks and navigated GCP's. These angles 
and the time correction are saved in an updated INF file W’hich will be used to navigate 
the image in the Forward process. 

5. Forward 

This routine produces a forward navigation of the image, as detailed in Appen¬ 
dix B. It utilizes the latest INF file and creates a N'AV (navigation) file. The navigation 
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Fig. 17. Twiddle procedure 


44 






file created contains tiie navigational values for every 16th line and pixel for the entire 
image. This file is cr*fltcu so ihat a latitude/longitude grid may be overlaid on a satellite 
image. Also this file provides a computationally quick method to obtain 
latitude/longitude values at different pixel locations when lower accuracy will suffice. 

6. Paint 

The Paint routine has two purposes. One is to validate the satellite image nav¬ 
igation and the other is to display sub-scenes of the satellite image. To aid in validation 
of the navigation, Paint performs a number of useful functions on full resolution sub¬ 
scenes of the image. Coastlines and political boundaries may be overlain on a sub-scene 
as well as latitude-longitude grids. The navigation values of any point can be displayed, 
allowing the operator to check how accurately the image has been navigated. Paint uses 
the INF, SCN, NAV and CAL files. 

7. Realmap 

This routine allows the user to map the sat lite image to one of 13 different 
projections. These projections are: 

Mercator 

North Polar Stereographic 

South Polar Stereographic 

Northern Hemisphere Lambert Conic Conformal 

Southern Hemisphere Lambert Conic Conformal 

Cylindrical Equidistant 

Azmuthal Equidistant 

.Modified Cylindrical uidistant 

Universal Transverse Mercator 

Transverse .Mercator 

North Orthographic 

South Orthographic 

Gnomic 

This enables the user to see the image in a convenient projection and utilize the data 
easily and effic..ntly. 

8. Utilities 

This routine enables the user to examine and manipulate the data or informa¬ 
tion files created throughout the Avian process. Landmarks can be deleted from the 
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LND flic if they are not located accurately enough. Information, such as the updated 
ascending node position and satellite attitude angles can be checked through this pro¬ 
cedure as well. 
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APPENDIX B. NPS SATELLITE IMAGE NAVIGATION 


A. COORDINATE CONVERSIONS 

Navigathg satellite images is equivalent to developing a mapping function from 
satellite image coordinates into geographic coordinates. This mapping is a complicated 
function and while it is possible to undertake in a single step, it is easier to understand 
with steps using intermediate coordinate systems. This outline of satellite image navi¬ 
gation is taken from Burks (1988). 

1. Coordinate sj’stems 

Developing the mapping from satellite image coordinates to geographic coordi¬ 
nates (and also the reverse process) involves the following six coordinate systems: 

1. satellite image coordinate system 

2. scanner coordinate system 

3. satellite orbit plane coordinate system 

4. geocentric orbit plane coordinate system 

5. geocentric coordinate system 

6. geodetic coordinate system 

For brevity the 'satellite orbit plane coordinate system' will be referred to as the 'satellite 
coordinate system' and the 'geocentric orbit plane coordinate system' will be referred to 
as the 'orbit plane cc ordinate system'. The first and the sixth coordinate systems in the 
list are two-dimensional coordinate systems. The three-dimensional coordinate will be 
in cartesian form, allowing easy matrix manipulations. Occasionally, spherical coordi¬ 
nates will be referenced and latitude, longitude and radius coordinates will be used. 

In the cartesian system, the equatorial (zero latitude) plane coincides with the 
cartesian x-y plane with the x-axis pointing through zero longitude, the y-axis at 90 de¬ 
grees from the x-axis, in the clockwise direction. The z-axis is normal to the equatorial 
plane, pointing in the direction of increasing longitude in a dextral sense, thereby defin¬ 
ing a right-hand coordinate system. 

a. The satellite image coordinate system 

This coordinate system is defined by the image scanning instrument. The 
origin (1,1) is the first pixel observed on the first line of the image. The first line of an 
image is arbitrar>, howe\er in Avian the first line is chosen to coincide with the earliest 
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time (fig. 18). Note that for a descending orbit, the first line will be shown as the 
northernmost point on the image; while for an ascending orbit, the first line will be 
shown as the southernmost point on the image. The x-coordinate is the pixel number, 
increasing along a single scan line. The y-coordinate is the scan line number, increasing 
along the image. As mentioned earlier, this system is a two-dimensional system with no 
spherical (or polar) coordinate counterpart. 

b. The scanner coordinate system 

This coordinate system is generally aligned with the satellite coordinate 
system, described later in this section. Its origin sits in the scanning instrument on its 
rotational axis. The x-axis is normal to the scanners' rotational axis pointing in the di¬ 
rection of the spacecraft motion. The z-axis points straight downward in the general 
direction of the center of the Earth. The y-axis completes a dextral (right handed) co¬ 
ordinate system. 

c. The satellite coordinate system 

This coordinate system has its origin at the instruments scan axis. The x- 
axis points in the coordinate of the motion of the satellite. The z-axis points through 
the Earth's center and the y-axis completes a dextral coordinate system. The x-z plane 
coincides with the instantaneous satellite orbit plane (fig 19). It differs from the scanner 
coordinate system by a rotation through the satellite attitude angles (roll, pitch and 
yaw). 

d. The orbit plane coordinate system 

This coordinate system is geocentric, with the x-axis pointing through the 
perigee of the orbit. The z-axi: is normal to this plane and points in the direction of the 
satellite's angular velocity. The y-axis completes a dextral coordinate system. In its 
spherical coordinate analogue, the longitude is called the true anomaly. The center of 
the Earth is one of the assumed elliptical orbit foci (fig. 20). 

e. The geocentric coordinate system 

This coordinate system is based on the Earth. The z-axis is in the direction 
of the Earth's rotational angular velocity. The x-y plane is normal to the z-axis, an¬ 
chored at the Earth's center. The x-axis points through zero longitude and the y-axis 
completes a dextral coordinate system (fig. 21). 

/. The geodetic coordinate system 

This coordinate system is the standard latitude, longitude system. It is 
based on the Earths' idealized oblate spherical surface. The third dimension (height 
above the Earth's surface) is not used for satellite image navigation (fig. 22). 
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Fig. 18. Chosen origin of AVHRR imagery in Avian (Rao et al., 1990) 

2. Geodetic to geocentric coordinates 

Geodetic latitude and longitude are the end product of the full satellite image 
navigation and typical of normal map presentations. However, it is connected with tlie 
Earth's surface. This conversion changes the origin of the coordinate system from the 
Earth's surface to the Earth's center and changes the reference suiface from an oblate 
spheroid to a sphere. This converts the geodetic latitude {(f)') and longitude (/T) to the 
geocentric latitude (d>) and longitude (.^) as shown in equation (B-1). 

<f) = atan{ tan </>' x (1 - J^)) 

;. = /l' (5-1)- 

Where f is the flattening of the Earth's surface (approximately lj291). The Earth's ra¬ 
dius (r) is a function of (p': 
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Fig. 19. Satellite coordinate system (Escobal, 1965) 
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1.0 



( sn\^<p) 

~V~' 


(5-2). 


where a is the semi-major axis and b is the semi-minor axis of the Earth's surface. 
Having the latitude, longitude and radius, a cartesian geocentric position vector for a 
point on the Earth's surface can be constructed. 

3. Geocentric to orbit plane coordinates 

This conversion is a well-defined three-dimensional rotation using the orbital 
elements at a given time. The rotation matrix G is given as: 


G = 


{cp X cn — sp y. sn y. cf) 

( —sp X c/i —cpx sn X Cl) 
{sn X si) 


{cp X sn + sp X cn X ci) 

( —sp X sn + cp X cn X cf) 
{—cnx sf) 


{sp X si) 
{cp X sf) 
ci 


(5-3) 


where: 

cp cosine of the argument of perigee 

sp sine of the argument of perigee 
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Fig. 20. Orbit piano coordinate system (Escobal, 1965) 


cn cosine of the geocentric ascending node longitude 

sn sine of the geocentric ascending node longitude 

ci cosine of the orbit inclination 

si sine of the orbit inclination 

This rotation is a combination of rotating through each angle individually, in the correct 
order. 

The reverse transformation multiplies an orbit plane coordinate position vector 
by the inverse of the matrix Cj. Since the rotation matrix is orthogonal, this is equal to 
multiplying by the transpose of the matrix. 

4. Orbit plane to satellite coordinates 

This tran.sformation is one of the more dilTicult used, as it involves both a ro¬ 
tation and a translation. Figure 23 depicts a satellite to geocentric transformation 
(Saufley, 1982). In this figure, c and d are angles which, if known, determine the position 
of a pixel on the ground. Both of these coordinate systems use two axes to define the 
orbit plane, differing as to \shich two they use. The position of the reference longitude 
in the orbit plane also differs by a rotation in the orbit plane through the argument of 
the perigee. 






Fig. 21. Geocf'jilric coordinate system (Escobal, 1965) 

Thus the full rotation is a multiplication of these two matrices (eqn. B-4). 


y?(y) = AT-(.v) 
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A' = 
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-1 


-1 
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0 


0 0 1 


thus, 
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—cw 
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(5-4) 


(5-5) 


(5-6) 


(5-7) 


52 








Fig. 22. Geodetic coordinate .system (Escobal, 1965) 
where 

cw cosine of the argument of perigee 

sw sine of the argument of perigee 

To translate the centers of the coordinate systems, the plane triangle defined by 
the center of the Earth, the satellite and the observed point is used. Three vectors from 
the sides of this triangle arc: 

S the satellite vector, which is the geocentric position of the satellite 

P the pixel vector, which is the geocentric position vector of the observed 

point 

V the view vector, which is the vector from the satellite to the observed point 

From these relationships you obtain: 

5 + P = r (Z? - 8) 

The pixel vector (P) and the satellite position vector (S) are known, so once the two 
vectors are rotated into the same coordinate system the third can be found (Fig. 24). 

The reverse transformation multiplies a satellite coordinate position vector by 
the inverse of the matrix (R). Since the rotation matrix is orthogonal, this is equivalent 
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PRIME 



Fig. 23. Satellite to geocentric transformation (Saufloy, 1982) 

to niultiphing by the transpose of the matrix R. The vector addition also still holds, 
except that the satellite and pixel vectors are known, giving: 

P=V-S {D-9) 

5. Satellite to scanner coordinates 

The scanner is solidly anchored onto the satellite, but the satellite may not be 
lined up in orbit. The small rotations from a true attitude are the satellite attitude angles 
(roll, pitch and yaw) which represent rotations around the x-, y- and z-axes respectively. 
Converting from satellite to scanner coordinates thus involves rotating the position 
vector through the roll angle, the pitch angle and then the yaw angle. Each rotation is 
represented by a rotation matrix R, P and Y respectively.. The total rotation is the 
multiplication of these matrices: 

A{x') = R X P X Y{x) {D-IO) 
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where 

cr 

cosine of the roll angle 

sr 

sine of the roll angle 

cp 

cosine of the pitch angle 

sp 

sine of the pitch angle 

cy 

cosine of the yaw angle 

sy 

sine of the yaw angle 


The reverse transformation multiplies a scanner coordinate position vector by 
the inverse of the matrix A. Since the rotation matrix is orthogonal, the inverse is equal 
to the transpose of the matrix A. 

6. Scanner coordinates to image coordinates 

This coordinate transformation is defined by the sensor model. The instrument 
scans the image and each pixel and line is mapped into a unit vector in the scanner co¬ 
ordinate system. If the scanner moves only in one plane the line number is a function 
of time. This is completely instrument dependent. 

The AVHRR radiometer scans within a single plane (Fig. 17), where the scan 
angle (scanL) is a linear function of the pixel number: 

scant = 9.4471 ne®"* x pixel no - 0.967856654 (5-15) 

where the zero angle points straight downward. Thus a unit vector in the scanner co¬ 
ordinate system can be constructed b>; 


.v = 0.0 


y = cos{scanL) 

z = s\niscanL) {B-16) 

These are the scanner coordinates of the pixel and the navigation proceeds from there. 

The scanner coordinates of the pixel are independent of the image line number. 
Since the scanner moves from line to line purely with the satellite motion, the line 
number is a strict function of time: 

^ ^0 + ^ — 1) (B—17) 
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where /o is the time for the first line in the image and the units are in seconds. Note that 
the AVHRR radiometer scans six lines per second. When navigating in a given direc¬ 
tion, the line number yields the time at which to calculate the orbital elements. 

B. ORBITAL ELEMENT MODEL 

An orbital element model predicts the orbital elements of a satellite at a certain time, 
given the orbital elements at an epoch time, usually once a day. The accuracy of this 
model is paramount in satellite image navigation, especially for polar-orbiting satellites 
which pass relatively close to the Earth at high rates of speed (making two revolutions 
per day), as opposed to geosynchronous satellites which orbit over one area or spot. 

The parameters used as orbital elements may vary according to the problem under 
investigation and the programmers discretion. A set always contains six parameters and 
any set of six can be converted to any other set of six parameters. An example of orbital 
elements is given in Figure 25 (Saufley, 1982). The six parameters used in Avian are: 

e the eccentricity is a measure of the difference of the satellites orbit from a 

circle 

i the inclination is the angle between the satellite orbital plane and the 

Earth s equatorial plane 

<y the argument of perigee is the angular measure of position of the orbit 

perigee (the point closest the center of the Earth) along the orbit relative 
to the Earth's equatorial plane 

V the true anomaly is the angular measure of the satellite's position relative 

to the orbit's perigee along the orbit 

a the semi-major axis length is half the distance along the long axis of the el¬ 

liptical orbit 

fl the longitude of the ascending node is the geocentric east iongituc> of the 

point where the orbit crosses the equator from south to north 

The angles measured are from the center of the Earth which sits at one of the foci of the 
elliptical orbit. The semi-major axis is measured from the center of the ellipse, not the 
center of the Earth. 

The orbital elements arc normally given at an epoch time which is not the time that 
is needed for the navigation, so the time will need to be updated. For Avian, the satellite 
ephemcris data is received with epoch times at OOOOZ each day. These Naval Space 
Surveillance Center elements include: 
eccentricity 
inclination angle 


^0 

4 
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T Epoch (date and time) 
a Semimajor Axis 
€ Eccentricity 
i Orbitol Indinolion 
(j) Argument of Perigee 

rj/ Mean Anomaly 

a Right Ascension of Ascending Nodp 



ASCENDING NODE 

A-Apogee 
P- Perigee 
S-Sotellite Position 
AN-Ascending Node 


Fig. 25. Orbital elements (Saufley, 1982) 

Wo argument of perigee 

Mo mean anomaly, wliich is the angular position of the satellite relative to the 

orbit's perigee, along a circle circumscribing the orbit 

/?o anomalistic mean motion (rad/herg) is the angular velocity of the satellite 

along the circumscribing circle 
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n' orbital decay rate (rad/her/herg) is not an orbital element, it measures the 

acceleration in the mean anomaly due to non-inertial and non-perturbation 
effects 

fto geocentric longitude of the ascenv>mg node 

Projecting these orbital elements forward in time is not difficult Four do not change 
with time: 


i=‘o 

<0 = (Oq 


a = 4~t^ (5-18) 

where the latter quantity is derived from Kepler's laws of orbital motion. The longitude 
of the ascending node only changes due to the Earth's rotation: 


(5-19) 

where is the Farth's rotation rate. The mean anomaly of the satellite is: 

il/=3/o-fV/ (5-20) 

where dt = T - Converting this to the true anomaly requires using elliptical ge¬ 
ometric steps. Firsi, find the eccentric anomaly (E) by solving Kepler's equation: 

M = E — e X siiiE (5 — 21) 

The solution is a zero order Bessel's function of the first kind. Following Smith (1980), 
expand the solution in terms of eccentricity a id truncate it at an approximate term. For 
this case, keep up to the third order. The result is: 

E = M+ sin(.\/ X e) 4 sin(;l/) cos(M x e^) + 

{sin.M - (-^) X sin(3 x A-/) x e^) + 0(e'’) (5 - 22) 

since e is quite small, the error is negligible. The true anomaly is found by: 
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sm V = ■ 


(^(1 - e^) X sin E) 
(1 — c X cos£) 


[ cos(£ — ff)] 
(1 — e X cosE) 


(5-23) 


Thus the six orbital elements have been obtained. 

There are other problems to be resolved. The Earth's gravitational potential is not 
spherical so the gravitational field exerts a force on the satellite. Since the Earth's shape 
is well known, these perturbation effects can be calculated. However, the results show 
that each orbital element depends on every other orbital element, thus creating a non¬ 
linear problem. Fortunately, if integrated over one orbit, the eccentricity, inclination 
angle and semi-major axis do not change, which makes the problem less complex. The 
orbital elements with the perturbative effects (to the first order) included are given by: 


(? = C0 


i = iO 


n = • 


(»o^ 


1 -i- 


(y)x J2X(1 -^Yt^x(l-(-J-)xsm 2 /) 


[«^x(l-eyj 


a = ti 



M = A/g 4- ndi 


OJ = COg + 


[(-|-)x72x(2-(y)x5/«2/)] 


[a^x(l-eY] 


dn 


( -^ X ^2 X cosi) 

£i.Qs--pi—--(B-24) 

la x(l-c) ] 


where J; is the first harmonic of the Earth's gravitational potential. Here, n is the mean 
motion constant, not the projected anomalistic mean motion. In the unperturbed case. 
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these two values are equal, eliminating any problems. In this case however, more care 
must be taken. The semi-major axis is derived from the mean motion constant, not the 
anomalistic mean motion. Deriving the true anomaly from the mean anomaly as de¬ 
scribed earlier, completes the orbital elements set. 

Another complication is resultant from other non-perturbative forces on the satel¬ 
lite. These forces include atmospheric drag, gravitational influence of other bodies and 
radiation pressure. All of the effects can be combined into the orbital decay rate given 
in the epoch orbital elements which results in a correction derived from observations. 
This factor only changes two equations, the mean anomaly and the semi-major axis; 


a = n 



(4 X a X «') 

-- X nxdt 


M = Mq ndt n'dt^ 


(5-25) 


Unfortunately, this introduces a non-linearity between the calculation of the mean mo¬ 
tion constant and the semi-major axis. The two quantities can be derived iteratively. 
First, calculate the semi-major axis using the equation: 

a = (5-26) 

Use this value to calculate the mean motion constant of equation (B-25). This can then 
be used to calculate the full semi-major axis length according to equation (B-25). Iterate 
over these two equations until the result converges to satisfactory limits. Twice is suf¬ 
ficient in this case. The rest of the orbital elements follow straightforwardly the 
equations of the perturbative case., 

C. REVERSE NAVIGATION 

The reverse navigation of an image converts the latitude and longitude of an ob¬ 
served point to its image line and pixel number. The process steps through the coordi¬ 
nate systems from the geodetic to the satellite image coordinate systems (in reverse order 
of the list in section A of this appendix). The reverse transformations are given, but 
there are several problem.s. The time at which the pixel was observed is unknown and 
the orbital elements at that time are unknown also. Fortunately, the orbital elements 
are veiy nearh constant with time or var} nearly linearly with time, so the solutions can 
be found iteratively. That is, assume a time, calculate the orbital elements, then update 
the time and run through the cycle again. 
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Examining the problem more closely, the unknown time means that the satellite 
position is unknown, which affects the conversion from plane coordinates to satellite 
coordinates. The rotation part is fine using assumed orbital elements, but one side of 
the triangle is unknown. The three sides are the view vector (V), the satellite vector (S) 
and the pixel vector (P). The three angles are the nadir angle (n), the zenith angle (z) 
and the scan angle (s). Of these six quantities, only the pixel vector is knowm. 

For a first approximation, assume that the true anomaly for the satellite is the same 
as the true anomaly for the pixel. This is a zero satellite attitude assumption, which is 
normally a good assumption as the satellite attitude normally varies only small amounts. 
Using this true anomaly (v sat), the satellite's orbit plane spherical coordinates are: 

longitude = vsat 


latitude = 0 


radius = ax 


(1 + c X cos{vsai)) 


(5-27) 


These can then be converted into cartesian coordinates, the satellite vector obtained, and 
from equation (B-8), the view vector is calculated. 

Continuing the process eventually yields a line and pixel number for the point, 
howe\er to iterate successfully, the time needs to be updated after each iteration. The 
satellite true anomaly (v sat) varies smoothly and monotonically with time, providing the 
needed informatic.:. Reversing the orbital element model yields: 


£ = atan 


( y '(1 — e*) X sin v) 
( cos V + e) 


M = E — esinE 


[B - 28) 


This mean anomaly is uncertain by a factor of 2n, that is, the number of satellite orbits 
since the epoch. Use the current estimate of the mean anomaly to get the integral mul¬ 
tiple of 271 required and add it to the calculated mean anomaly. The time since epoch 
is derived by solving the quadratic equation: 


M = .Vq -t- ndt -f fi'dt^ 


{B - 29) 
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a = n' [B - 34) 

the time can is then obtained from these equations through iteration. 

The true anomaly currently known is for the pixel. Accuracy requires an estimate 
for the true anomaly of the satellite. With an estimated pixel number, the offset between 
the satellite and pixel true anomaly can be calculated quite precisely. 

Examine Figure 24, the plane triangle of a satellite's view. The scan angle (s) is the 
arc between the sub-satellite point and the pixel. Only the satellite attitude angles keep 
the orbit plane from being normal to the triangle. The arc normal to the orbit plane 
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from the pixel to the orbit track completes a spherical surface triangle. In orbit plane 
coordinates, the arc length of the normal arc (b) is the orbit plane latitude of the pixel 
which is known. The offset between the pixel and satellite true anomalies is the arc dv. 
If the attitude angles are zero, the offset is zero and the satellite true anomaly is equal 
to the pixel true anomaly. 

To solve for the offsei, one more quantity needs to be known about the spherical 
surface triangle. The angle y is easy to find, as it is related to the longitude of the pixel 
in the satellite coordinate system. 


y = 


TT 

0 


— aian 


xsat 

ysai 


[B - 35) 


The X and y satellite coordinates can be derived by taking the first few steps of the for¬ 
ward navigation algorithm (Section D), converting the pixel number to satellite coordi¬ 
nates. With y known, the rules for a right spherical surface triangle give: 


dv = asin 


{tan b) 
(tail y) 


{B - 36) 


However, the orbit plane latitude (b) and y are both signed values, so there are four cases 
to be considered: 


Table 4. SIGN OF DV DEPENDING ON LATITUDE AND GAMMA 


Latitude 

r 

sign of dv 

pos 

<pi'2 

pos 

pos 

> pi 2 

neg 

neg 

<pi'2 

neg 

neg 

mmmBmm 

pos 


Thus the satellite's true anomaly is: 

vsat = vpix -f dv {B — 37) 

With this final equation, the reverse navigation is well-defined. 
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D. FORWARD NAVIGATION 

The forward navigation of an image converts the image line and pixel number of the 
observed pixel to its corresponding latitude and longitude. The proces. goes through the 
coordinate transformations from the satellite image coordinate system to the geodetic 
coordinate system, as shown previously. 

Just as with the reverse navigation, there is missing information that needs to be 
known before the navigation can be accomplished. The position of the pixel on the 
Earth's surface is unknown. 

The key to obtaining this unknown information is in the triangle formed by the 
center of the Earth, the satellite and the pixel. The satellite vector (S) from the Earth's 
center to the satellite is known. The direction of the view vector (V) from the satellite 
to the pixel is known. The pixel vector (P) from the Earth's center to the pixel is un¬ 
known. The length of the view vector can be calculated if the Earth's radius at the pixel 
is known. The nadir angle (n) between the view and satellite vectors is the arccosine of 
the dot product of those two known vectors: 

r {.xvicw X xsa! + wiew xysai + zview X zsat) "1 -.ns 

n = acos\ -=- - - {B — 38) 

where r is the length of S (i.e., the satellite's geocentric radius). From the law of sines: 


a _r__ 

(sin n) (sin 2 ) 


{B - 39) 


where a is the Eart; radius and z is the zenith angle plus ninety degrees, and 


V _ r 

(sins) (sine) 


{B - 40) 


where v is the length of the view vector and s is the scan angle. The sum of s, n and v 
is equal to n. Combining this with equations B-39 and B-40 gives the length of the view 
vector V. With the view vector known, the pixel vector is given by: 


P = 5+ K 


(B-41) 


The Earth radius at the pixel is still unknown. Since the Earth's radius is constant within 
25 km, its effect is second order and an iteration can solve the problem. This is accom¬ 
plished by assuming any initial latitude for an approximation of the Earth's radius. 
Then iterate over the full conversion, using that latitude. When the pixel position is 
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calculated, update the estimated latitude and iterate again, continuing until the needed 
accuracy is obtained. Then the forward navigation process can be completed. 
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appendix c. code for reading vvvs 


Program WVS 

c* 

C* Purpose! 

c* Tills purpose of this routine Is to let the user select a section 

c* of VVS coastline (by latitude range). The position of the coastline 

C* Is then written to a file (lat lon.dat) for later use or plotting. 

C* The density of the coastline (T.e., every point or every other point, 

C* etc.) can be chosen Is the subroutine show_seg using the skip/ / 

C* command. 

c* Presently, this routine Is In developmental form and changes 

c* are imminent. This program worked well for the west coast of the 

c* United States. There were problems outside of the U.S. boundaries 

C* possibly due to Internalonal boundaries. This problem will be 

C* dealt with at a later time, as this program Is in its early 

C* development and for our purposes, only the U.S. west coast was 

C* necessary. 

C* 

C* Variable naming philosphy: 

C* For the most parts the same name Is given to variables In this 

C* computer code as used In VVS specifications. 

C* 

c* Major variable list: 


c* clat . Floating point equivalent of latcell. 

c* cion . Floating point equivalent of loncell. 

c* Irec .. Logical record i (48 bytes). 

c* ipos . Used to store location of byte pointer and 

c* pass this value to other subroutines. 

c* nbytlft . Number of bytes left In present physical record 

c* to where the next cell (C) starts. 

c* nrecleft . Number of logical records until next cell (C). 

e* npiectosklp . Number of physical records to skip before getting 

c* to next cell. 

c* prec . Physical racotd H (9600 bytes). 

c* 

c* Problems: 


c* The routine as been tested on any general cases. Only one data tape 

c* "as utilized. The data tested didn't have any extra attributes, etc. 

c* In addition, the bad programming technique of common blocks was 

c* utilized for my convenience, 

c* 

Implicit none 
integer*4 blksize 
parameter (blksize - 9600) 

lnteger*4 cellnum, Irec, i,j,channel, Istat.slzeread 

lnteger*4 datadec, Ibeg, latcell, Ingcell, Ipos 

lnteger*4 clat,cing, cellstart 
integcr*4 minlat.maxlat 

lnteger*4 nfea rcc, ncellhead, nrec left, nseg rec 

lnteger*4 nleft, : nprectoskip, ntext, origdec 

lnteger*4 nfeaincell, nsegincell 

lnteger*4 bytepos, size, status 

logical tape 

external bytepos 

character*9600 chrbuf 

character ans*l, ass*10,format*ll,ccllflag*l 
byte buf(9600) 
logical echo, dataflag 
equivalence (chrbuf,buf) 

Integer prec 

common /physrec/chrbuf, prec 
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data ptec/O/i nle£t/0/ 
data delaflag/.true./ 

.. define functions 
integer**! skip_precs 

print '(a,S)*i' Read data from tape (t) or disk (d) ->' 

read '(a)',8ns „ 

if t»nv eo. 't' .or. ans .eq. T )tnen 

J?^;f:;“r?l?Tchl*:nind(cHannef. istat, 1 ^SS^nr^^pt'dlHe 
if (istat .ne. l)then 

print open status Istat 

stop 'Problem opening tape drlvet aborting 

print *)' Routine opened taped channel successfully 
end If 

*^*open (unlt.lO.flle.'dfskSlltiscratch.motelllwvs.dBf, 

\ ACCESS*'DIRECT',£orm»'UNrORHATTED',recl-blkslzc/4, 

2 iostBt*status,status-'OtO') 

l£*(status'°ne. Olprlnt Status problem opening diskfile' 
endl f 

print '(a)',' User Has option to see every bit of data. But' 

Lint 'L)'.' 1 vould advise you to Just say HO to check data 
print '(8,$)',' Check data (y/n) •>' 

read '(a)',ans , .vm „ 

if (ans .eq. 'Y' -or. ans .eq. 'y')then 
print *,' enter cell number starts’ 
echo • .true, 
read *,cellstatt 
else , 

echo • .false, 
endi f 

print *.' Echo values' 
print ‘teelio 

print ' (a,/|a,/.a,S)', r 

1 ' Enter latitudes rangefpositlve for North) , 

2 ' For example, •>32, 38', 

3 ' Enter latitudes ->' 

read *,rolnlat,maxlat , , . i i , 

print *,' Using following lat,lon boundss ,minlat,maxl8t 

open (unlt.l2,fl!e»''l8t Ion',Bccess.'sequential', 

( form.'formatted',status-’new') 

open (uni t.l5,file-'SH0RE.DAT',access.'sequential , 
i form.'formatted',status-’new') 

print '(a)’,' Opening file <lat_lon.DAT> for data output 
V , 

... read at least the first four file header logical records 

print reading first physical record' 
call get_data(channel, echo) 


.... decode first file header record 
print 700 ,chrbuf(1sA8) 
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read (chrbuf(1i48),1000)origdec,datadec 

1000 form8t(38x,J5,i5) 

piint origdec »',orlgdec 

print datadec a'ldatadec 

c 

c.... decode second file header 
c 

print 700,chrbu£(49!97) 
c 

c.... decode third file header 
c 

print 700,chrbuf<98!l45) 
c 

c.... decode fourth file header 
c 

print 700,chrbuf(146:193) 
read (clirbuf( 146:193), 1001 )nlext 

1001 format(41x,14) 

print ntext -'.ntext 
print 705,ntext 

705 formate Skipping <',15,'> text records') 
c 

1 . 193 
C 

C.... Read Another piiyslcal record from tape: 

C 

5 continue 

If (prec .gt. l)then 
print 830,pree*l 

830 formate reading physical record *',15) 
call get datafchannel, echo) 
endlf 
c 

c.... process data 
c 

10 continue I Increment pointer to data within physical record 

c 

c.... Perform data search: first check to see If lat h Ion are In bounds, 
c.... If not jump to the next cell. Use Information In cell header to 
c.... determine how many logical records roust be skipped before next cell, 
c 

715 format(lx,a) 
read(chfbuf(l:l),I025)cell£lag 
do while (cellflag .ne. 'C') 

print 7i6,chrbuf(1:1447) 

716 formate cell mismatch:',a) 

1 » byteposfchannel,1,1,echo) 
re3d(chrbuf(l:i),1025)cellflag 
end do 

1025 format(al) 

c 

leadfchrbuf(1:1447),1010)celIflag,ncellhead,cellnum, 

4 lngcell,latcell, 

4 nfealncell,nsegincell,nfe3 rec,nseg rec 

iOlO format(al,lx,11,16,18,17,15,15,17,17) 

706 £ormat(6x,a) 

If (cellflag .eq. 'C' .and. nfearec .ne. 0)then 
clng « Ingcell/origdec 
clat • latcell/orlgdec 

If (rlat .ge. roinlat .and. clng .le. maxlat)then 
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print 715ichrbu£(l:l+47) 

call show_le8(c«llnum, channel,!, clat, clng, datadec, 

& “ echo, nfeaincell, nseglncell) 

print 715,chrfauf(l!l+47) 

1! (1 .gt. 9600)then 
1 . 1 

goto 5 I read new physical record 

end! f 

goto 10 I process next logical record 

else 
c 

c.... Determine II logical records to skip to next 

C-... Must alvays move byte pointer "1" at least oi.e logical record. 

c 

nrec left ■ ncellheod*nfc8 rec«nseg rec 
c 

C-... Skip unwanted logical records by Incrementing pointer (1) 
c 

1 • byte posfchannel, 1, nrec left*1,echo) 

1F(I.E0.0)TIIEN 

PRINT I .0 HERE' 

GOTO 5 
end! f 
goto 10 
endl f 
end! f 

I • bytepos(channel,1,1.echo) 1 point to next logical record 

goto 10 I process next logical record 

c 

700 format(lx,a) 
end 
c 
c 
c 

subroutine get data(channel, echo) 
c* 

c* Purpose: 

c* Routine reads one physical record from tape or disk. Assumes 

c* that logical device unit or channel • 10, when Its from disk, 

c* If an installation every has a tape unit device number of "10", 

c* this routine would become confused, 

c* 

c* Argument list: 

c* Istat . Status of tape read (1.success). 

c* 

c* Problems: Uses common block which Is a bad programming practice, 
c* 

lnteger*4 blksize 
parameter (blksize • 9600) 
lnteger*4 channel, Istat, sireread 
ch3racter*9600 chrbuf 
byte buf(9600) 
equivalence (chrbuf,buf) 

Integer prec 
logical echo 

common /physrec/chtbuf, prec 
c 

prec « prec*l 

If (mod(prec,10).eq.0)prlnt Physical record *',prec 
If (channel .ne. 10) then 
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prim Reading Iron tape ' 

call mtread(ctiannol,buC,blkslze,lstatislzeread) 
if (slzecead .ne. biksize)then 
print Number of bytes read a',sizeread 

print *i' Should have read a'lblksize 
print 7 what''s up 7' 

prec « prvc-1 
endi f 
else 

read (channel,RECaprec,ERR.99)buf 
end! f 


if (echo)then 
do 1.1,9600,^6 

print 716,1,chtbuf(1:t447) 
write <15,716)i,chrbut(i:lei?) 
end do 
end! f 

716 formate G.D. ',14,lx,a) 


c* 

c* 

c* 

C' 

c* 

c* 

c* 

c* 

e* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c» 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c8 

c* 

c* 


return 
99 print * 
return 
end 


error reading recordsprec,channel 


subroutine show_fea(cellmin, channel,!, clat, clng, datadec, 

& ” echo, nfeaincell, nsegincell) 

Purposes This subroutine is used to decode Inforsnatlon from the feature 

record. After needed information is decoded it is viltten to a file, 
via other subroutines. 

Major Variabless < 

contl_arr . This array stores the edge code (see pp. 14-15 of 

VVS) as a function of "feanum". 

feaflag . Character variable that should be 'FEA' when 

decoded. 

feanum . Don't understand what this is yet. 

nfdatarec .The I of feature data records that follow the 

feature header, usually 1. 

segdlr_arr . This array stores the direction of the segment 

(either "F’’,"R" or "D"). The direction is stosed 
as a function of segment ID. 

seg_assoc . This array stores the feature number "associated" 

with a given segment. Thus, the indice of "seg_8SSoc" 
is the segment number, and the value is the 
associated "feature number". 

Algorlthm; 

1) The cell header gives the total number of feature records. The 
first step, therefore, is to process all feature records. 

2) Vhen processing feature records, store in array the "edge code” 
associated with a given feature number. 

2) Next, the number of segment records is used as a loop limit. That 
is, process each group of segment records (and associated vertices) 
until no segment records are left. Note that there can be more 
segment records than feature records. 
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c* 

c* 

c* 

c* 

c* 

c* 


Notes- Msnual overltjo. H has been found that occasional the J’*® 

Notes. t,, nfdatatec or Number of Feature DATA RECotds. 

Thus, the routine skips to where It thinks the next feature o 

segment record should reside and finds a different 

expected. This routine allows one to manual overlde this proble . 

implicit none 
charactpr*9600 chrbuf 
integer prec 

common /physrec/chrbuf, prec 

Integer*^ cellnum, clat.clng, channel, datadec 
Integer** 1 
logical echo 

... local variables . . 

Integer segnum_arr(300), Index, 

Integer coiitl ,”cont2, contl_art(300), cont2_Brr(300) 
character*! segdlr 8rr(300)T seg_clir(6) 
character feaflag*3, featype*l,facs*5,segdir*l 
Integer** fealeft, by tepos,nfealncell,feanum.nseglntea 

Integer** nattr, nfdatarec 

Integer** nseglncell, nsdatarec, loop, segnum 
Integer** k,kstatt,j 


L.. Increment byte pointer to next record which should be a feature record 
c 

if (echo)then . , , . n 

print *,' 1 ,nfealncell,nseglncell,ptec',1,nfeaincell, 

(, nseglncell,prec 

print 715,chrbuf(l!i*A7) 
end 1 f 
c 

tea left • nfealrcell 
1 .“byteposfchannel,1,1.echo) 
kslart . 1 

do while (fea left .gt. 0) 
fea left - lea left - 1 
If Techo) print 715,chrbuf(i:1**7) 

read(chrbuf(l!l*47),1000)£cafl8g,£eanum,featype, 

g centl,cont 2 ,nsegInfea,nattr,nfdatarec 

c 

c... Store values need for later in arrays: 
c 

contl_arr(feanun.) • contl 
conl2 acr(£eanum) . contZ 
1000 form3T(83,17,al,5x,18x,ll,il,*x,i3,i3,i2) 

if (echo)then 

print 800 ,feaflag.feanum.featype,centl.contz, 

(, nseginfea,nattr,nfdatarec 

800 formatf' feaflag-',a3,' fcanum.',12,' featype-',al, 
h ' contl-',11,' cont2-',il,' nseginfea- ,1*, 

g ' nattr .',12,' nfdatarec.',13) 

end! f 
c 

If (feaflag .ne. 'FEA')then 

print *,' Next record Is not a feature record 
print 715,chrbu£(i:l-t47) 
stop ' Program aborting III' 
end! f 
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c 

c.... At this point the previous record most of been a feature record 
c.... Now skip past exti’ attribute records: 
c 

1 - kytepos(channel,i,nattr4l,echo) i skip feature data record, 
c 

c.... Decode feature data records ;nd store direction and number in arrays, 
c 


do j>l,nfdatarec I loop fot each logical record 

readichrbuf(1:1447),1050)(scg_int(k),k>l,6) 

1050 format(6(i7,lx)) ~ 

read(chrbu£(1:1447),1060)(seg_chr(k),k.l,6) 

1060 fotmat(6(7x,al)) ~ 

c 

c.... store segmr"t Identification number only for actual numbers: 
c 

do k«l,6 

1f<seg_lnt(k) .ne. 0)then 
kstaTt « kst3rt4l 
index - seg int(k) 
segoir arr(lnoex) • seg chr(k) 
end if 
enddo 
c 

1 « bytcposCchannel,1,1.echo) 
if (kstart .ge. 194)stop ' kstart is too large ' 
end do 
c 

c.... point to first segment record: 
c 

5 continue 

If (echo) print 715,chrbuf(1:1447) ! ECHO 

715 £ormat(lx,a) 
end do 


c 


c* 

c* 

c* 


c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


... The following records are SCCment records: 

... Pointer should be positioned at segment header record. 

do loop=l,nsegincell 

1f(echo)ptint loop»,nseginc€ll‘,loop.nsegincell 
call show_seg(cellnum, channel, clat,clng,contl_arr,conl2_arr, 
6 datadec.eclio, 1,loop,segdir_arr) ~ "" 

end do 
retuin 
end 


subroutine shov_seg(celInum, chn.clat,clng,contl err.contZ arr, 
h datadec.echo, 1 ,scgld,Sv>gdir_irr) 


Aigori thm;. 

1) Set value of seg verts • 0, this value will then be 
incremented until termination of routine or until seg verts 
is greater than some pre.set number of the lat.lon data arrays 
(l.e., usually 5000). 

3) S'art a loop to go through every segment record. 

2) Uoe parameter "skip" to determine the number seg verts 
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c* that vill be written to disk. This parameter is used since the 

c* resolution ot the world vector shoreline is much greater than 

c* needed. By modifying "skip", the user can skip vertices, thus 

c* decreasing the VVS resolution. 

c* X) "move_ptr" is used to move pointer past the blank segment vertices 
c* which'appear at the end of a group of segment records, 

c* 

c* Problems: Routine only works properly when skip^i, this is because 
c^ nrecs left parameter is not decremented properly otherwise, 

c* ~ 

implicit none 
character*9600 chrbuf 
integer cellnum, prec 
common /physrec/chrbuf, p'ec 

integer** clat,clnK,contl_arr(*),cont2 Brr(*),segid 
character*! segdlr_arr(*)” ~ 

Integer bytepos 

Integer chn,datadec,ptr, segnum, nverts, nsdatarec, nxfearecs 
c 

c.... Define local vatiables: 
c 

character segflag*! 
integer** feanun, Iflag 

integer** l,j,n,nrecs left, nverts left, seg verts,ians 
integer** counts, Istat, move_ptt,~skip, xseconds,yseconds 
real** invdiv, latarr(30000),”lon8rr(30000) 
logical eclio, bcg^seg 
data seg_vet ts/0/7 sk.ip/2/ 
save seg”verts 
c 
c 

beg_seg « .true, 
counts - 0 

if (eciio) print 715,chrbuf(1:i**7) 

715 format(lx,a) 

read(chrbul(i:l*47),1030)scgflag,segnum,nverts,nxfearecs, 

(< nsdatarec, feanum 

1030 format(a3,17,15,2x,i2,i5,i7) 
c 

c.... According to previous assumptions, segid =. segnum 
if (segid .ne. segnum)then 

print *,' segid '.segid,' does not equal segnum '.segnum 
return 
end if 

.... IF THERE ARE EXTRA SEOHENl RECORDS SKIP THEM: 

i - byteposCchn, 1, l*nxfearecs,echo) 
c 

c.... Pointer (plr) should be positioned at first verlicc pair: 
c 

nrecs_left « nsdatarec 
invdiv . l./(datadec*3600.) 
c 

if (echo) print 715,chrbuf(i!l+47) 
do wliile(ntecs left .gt. 0) 
c “ 

c.... Calculate number of vertice records left: 
c 

nrecs left « nrecs left - 1 
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nverts_le£t « nverts-counts 
.... Determine how many vertices to process: 


c 



1000 


1 

2 

3 


n » 4 I maximum vertices per record 

If (nverts left .It. 4)then 
n • nverts_left 

print nverts left*, counts',nverts left,counts 
endif 

move_ptr . (4 - n)«12 

do j*l,n I process four or less vertlces/rccord 

counts • counts + I 
if (modfcounts,skip) .eq. 0 )then 
scg_verts • scg_verts4l 

l£“(seg_vetts Tgt. 10000)stop ' seg_verts • 10000' 

If (seg'verts .gt. 7500)stop ' normal termination ' 
call closeup ( counts, seg verts, 

latarr(l), lonarr(l), Iatarr(5000), 
lonarr(5000) 

end! f 

print 715,chrbuf(l:l4ll) 

read (chrbuf(l:l4ll), 1000 )xscconds,yseconds 
format(i6,16) 

latarrfseg verts) • clat ♦ flo 8 t(yseconds)*lnvdiv 
lonarr(seg~verts) - clng 4 £loat(xseconds)*invdlv 
l£(beg segTthen 
Iflag * 0 
else 

iflag » I 
endif 
print *, 1 , 

seg_verts,latarr(seg vects),lon 3 rr(seg_verts), 
iflag, 

con tI_arr(feanum),cont 2 _arr(feanum) 


vritc(12,700) 

1 cellnum, seg_verts,latart(seg_verts), 

2 lonarrfscg verts), Iflag, 

3 contl arrfleanum), cont2 arrffeanum) 

700 format(2x,16,lx,15,f10.5,lx,flO.5,lx,3(11,lx)) 

beg seg * .false, 
endif” 

1 * 1*12 

If (1 .ge. 9600)tlien 
call get data(clm, echo) 

1.1 " 
end If 
end do 

1 * 1 4 move ptr 

if (1 .eq. 9600)print *..••...** 1 *',i 
end do 
return 
end 


c* 


c* 

Integer function skip prccs(channel,nprectoskip,echo) 
c 

Implicit none 
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Integer prec 
character*9600 chtbuf 
common /physrec/ chrbufi prec 
integer*^ nprectoskip,j,Istst.channel,nleft 
INTEGER I 
logical echo 
c 

do 30 j'l.nprectoskip 

call get data(channel,echo) 

30 continue 

skip precs • prec 
return 
end 
c* 
c* 
c* 

Integer function byteposfchannel, pointer, Irec, echo) 
c 

c* Purpose: 

c* This routine Is used to move the "byte pointer" to where 

c* the next logical record needed Is located. Dy setting the logical 

c* record size to zetoe, this routine can be used to place pointer In 

c* in next logical record. Vhen Irec Is greater than zeroe this routine 

c* will move the pointer to elthet the correct location In the present 

c* physical record, or the needed location In the next physical record 

c* according to where the next needed data should be located, 

c* 

Implicit none 
character*9600 chrbuf 
Integer prec 
logical echo 

common /physrec/ chrbuf, prec 

integer*4 channel, pointer, Irec, prectosklp, bytes_left 
lnteger*4 Istat,skip_prees, 1 
external skip precs 
c 

bytes_lcft • lrec*48 » pointer 
If (bytes_left .It. 9600)lhen 
bytepos'. bytes_left 
return 
else 
c 

c.... check to see If whole physical records need to be skipped 
c 

bytes_left » bytes left - 9600 1 I bytes in present record 

prectosklp = bytes left/9600 I I physical records to skip 

if(prectoskip .gt.~0)then 

prec • skip_precs(channel,prectosklp,echo) 
c 

c.... set pointer to new value in new physical record 
c 

bytes left • bytes left - prc.toskip*9600 
end If *" ~ 

endif 
c 

call get_data(ch3nnel, echo) 
bytepos » bytes_left 
c 

return 

end 
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c* 

c* 

c* 

subroutine closeup(prec,n,latl,lonl,lat2,lon2) 

Integer prec,n 
real latl,lonl,lat2,lon2 
c 

print Program processed:' vertices at prec: ',prec 

print lat,lon range from:',iatl,lonl,' to ',lat2,lon2 

close (15) 
close (12) 
close (10) 

stop ' normal termination ' 
end 
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